{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import glob\n",
    "import sys\n",
    "\n",
    "from scipy.interpolate import interp2d\n",
    "\n",
    "# Our modules\n",
    "sys.path.append('/home/jovyan/segtrax/source')\n",
    "from utilities import transform_coord\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proof of concept for vectorized indexing in xarray\n",
    "Create a DataArray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray (x: 3, y: 4)>\n",
       "array([[ 0,  1,  2,  3],\n",
       "       [ 4,  5,  6,  7],\n",
       "       [ 8,  9, 10, 11]])\n",
       "Coordinates:\n",
       "  * x        (x) int64 0 1 2\n",
       "  * y        (y) <U1 'a' 'b' 'c' 'd'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "da = xr.DataArray(np.arange(12).reshape((3, 4)), dims=['x', 'y'], coords={'x': [0, 1, 2], 'y': ['a', 'b', 'c', 'd']})\n",
    "da"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate index DataArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind_x = xr.DataArray([0, 1], dims=['i'])\n",
    "ind_y = xr.DataArray([0, 2], dims=['i'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray (i: 2)>\n",
       "array([0, 6])\n",
       "Coordinates:\n",
       "    x        (i) int64 0 1\n",
       "    y        (i) <U1 'a' 'c'\n",
       "Dimensions without coordinates: i"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "da[ind_x, ind_y]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Try it with ice motion data\n",
    "\n",
    "Read the ice motion data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:                 (time: 70, x: 361, y: 361)\n",
       "Coordinates:\n",
       "  * x                       (x) float64 -4.512e+06 -4.487e+06 ... 4.512e+06\n",
       "  * y                       (y) float64 -4.512e+06 -4.487e+06 ... 4.512e+06\n",
       "  * time                    (time) object 2019-01-01 00:00:00 ... 2018-12-24 12:00:00\n",
       "Data variables:\n",
       "    crs                     int32 -2147483647\n",
       "    latitude                (y, x) float32 29.896942 30.080757 ... 29.896942\n",
       "    longitude               (y, x) float32 -45.0 -44.8404 ... 135.15959 135.0\n",
       "    u                       (time, y, x) float32 dask.array<shape=(70, 361, 361), chunksize=(18, 361, 361)>\n",
       "    v                       (time, y, x) float32 dask.array<shape=(70, 361, 361), chunksize=(18, 361, 361)>\n",
       "    number_of_observations  (time, y, x) float32 dask.array<shape=(70, 361, 361), chunksize=(18, 361, 361)>\n",
       "Attributes:\n",
       "    version:       QuickLook data based on Ice Motion version 4.1\n",
       "    release_date:  May 2019\n",
       "    Conventions:   CF-1.4\n",
       "    citation:      Tschudi, M. A., Meier, W. N., and Stewart, J. S.: An enhan...\n",
       "    dataset_doi:   10.5067/O0XI8PPYEZJ6"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "motion_files = glob.glob('/home/jovyan/segtrax/data/icemotion*.nc')\n",
    "ds = xr.open_mfdataset(motion_files, concat_dim='time', data_vars='different')\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot one week of data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7fa43d3bdef0>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEWCAYAAAA3h9P4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOy9ebwlV1X3/V17V53hDn17SpM5nUCYDUEioMhDEAyoCI+CPhGFqMgkvPg++KqAPoDggD4iisoQTSBMxoj6gjwMRjCCL4gJYRZiQiaaJHR3uvv2nc45VbXX+8feVafu7Tuc7r7d95zb+/v51KfPqaq9a9fpe+p31tprryWqSiQSiUQiw4rZ6AFEIpFIJLIaUagikUgkMtREoYpEIpHIUBOFKhKJRCJDTRSqSCQSiQw1UagikUgkMtREoYosi4icKyKzImI3eiyRSOTUJgpVBAARuVNEnlq+V9W7VXVCVYuNHNdKiEhTRK4SkbtEZEZEvigiP7LknKeIyDdFZF5E/kVEzqsde3LYNy0idy7T/8Ui8plwfI+IvHaN8bxRRL4qIrmIvH7JsTNE5MMico+IqIjsHuD+nhvubU5E/l8R2V479m4R6YUfEuW27A8KEWmIyAfD/6+KyKVLjv+aiHwtfIZ3iMivrTGu3eFzmw+f7VOXHF9x3Mv01RSRq0XksIjcJyKvXHL8YhH5QrjWF0Tk4tXGFtm8RKGKjCoJ8G3gScAU8L+A60oREJGdwN+H/duBm4C/qbWfA64GVnowfwD4dGj7JOClIvLMVcZzG/DrwP9Z5pgDPg48e+3bAhF5BPBO4HnAA4B54G1LTvvD8ENiYoAfFP8G/Bxw33KXA54PbAOeDrxcRC5fpa+/Br4I7AB+E/igiJx2FOOu83rgQuA84MnAr4vI00NfDeBDwPvC2K4BPhT2R041VDVup/gGvBf/MF0AZvEP3N2AAkk45wbgd4DPhnP+Ef+wej9wGLgR2F3r86HA9cAB4Bbgp0/CfXwFeHZ4/SLgs7Vj4+H+HrqkzVOBO5fpax54eO393wKvHmAM7wNev8KxJHymu9fo4/eAD9TePxDoAZPh/buB3zmGz2cPcOka57wV+LMVjj0Y6JbjCPs+A7xkkHEv0993gMtq798IXBteXxaOS+343cDTT+Z3I27DsUWLKoKqPg//EPhx9b/O/3CFUy/H/1o+C/8Q+hzwLrzV8Q3gdQAiMo4XqQ8Au4CfAd4WfnEfgYi8TUQOrbB9ZZB7EJEH4B+kXw+7HgF8uXaPc8C3wv5B+BPg+SKSishDgO8H/nnAtsfL0rF/C//Af3DtnF8WkQPBJTaQpbYWIiLAE+l/hojIR0TkVbVx3a6qM7VmX6b/ma46bhF5lYh8JLzeBpxZP3+Zvr6iqvUcb19h8P+/yCYiClXkaHiXqn5LVaeBjwHfUtV/VtUcb3E8Opz3DLyV8i5VzVX1ZuDvgOcs16mq/rKqbl1hu2itQYlIirfsrlHVb4bdE8D0klOngckB7/UjYbwLwDeBq1T1xgHbHi9rjf2teJfZLrxr890i8oR1uO7r8c+Ed5U7VPUZqvqmAce16nFVfZOqPqN2LkvOH7ivyKlFFKrI0fDd2uuFZd6XD5/zgMfVLSPgZ4HT13tAImLwrsse8PLaoVlgy5LTtwAzrEEIAPg48AagBZwDPE1Efjkc/3otiOGJxzn+J9b6Ki2ZVceuqjer6v3hR8BH8SL9k8c5jpfj56p+TFW7K5y21md6NJ/5bO348fYV2eREoYqUrGca/W8D/7rEMppQ1Zcud7KIvGNJBNvsMg/v5doJcBV+4v7ZqprVDn8deFTt3HG8u3LF/mpcABSq+p4gBnuAa4EfBVDVR2g/iOEzA/S3Iqr6mVpfpVtr6dgvAJrAf63UDT4o4pgQkV8EXgU8JdzrSnwduEBE6lbNo+h/pgOPW1UPAvfWz1+mr4vC/3HJRQz2/xfZZEShipR8F/+AXg8+AjxYRJ4X5nhSEfk+EXnYcier6kt0cQTbxDIP7+V4O/Aw/NzawpJj/wA8UkSeLSIt4LX4OY9vgrfEwv7Uv5VWLaLsv8K+54bzTgf+B4vnUxYR7rGF/04loT9bO97CP7QBmuH9Srwf+PFgbY3jLbu/L+eGROQ5IjIRxnYZPqLvw6uMrX69RhibhGM/iw+C+GFVvX2VMaGq/wV8CXhd6OMn8OLxd4OMexneA/yWiGwTkYcCL8QHioAP3imAV4Txl9byp1YbY2STstHRHHEbjg14Fj6g4hDw/7B81N8v1c7/HeDdtfdPBW6rvX8IPlR7H3A//gFz8TqO97wwvg7eTVRuP7tkTN/EuyVvYHFU4qWhfX27oXb8h/CRjNP4sO6/BMZWGc+7l+nv52vHlx7TNe7vueH/Yw4fpr29duwzYVyH8eJ5+Rp93bnM9XeHY3cA2ZLP8B21th8DXlN7vzt8lgv4aM6nHsW4XwN8rPa+iV8icBj/Q+mVS/p6NPCFcK2bgUdv9PckbhuzSfiDiEQikUhkKImuv0gkEokMNVGoIpFIJDLURKGKRCKRyFAThSoSiUQiQ02y0QPYDOzcuVN379690cOIRCIjwBe+8IX9qnra8fRxjrS1gxvo3P30PqGqTz+e6200UajWgd27d3PTTTdt9DAikcgIICJ3HW8fHRzP5oyBzn0nd+083uttNFGoIpFIZMQQwA6ai2QTrECKQhWJRCIjhgANM6BSDWXp06MjClUkEomMGN6iOub0jiNHFKpIJBIZNeQoXH+bgChUkUgkMmJEiyoSiUQiQ81RBVNsAqJQRSKRyMgh0aKKRCKRyPAiQHoKCdWGp1ASESsiXxSRj4T320XkehG5Nfy7rXbuq0XkNhG5RUSeVtv/GBH5ajj21lpRuKaI/E3Y/3kR2V1rc0W4xq0ickVt//nh3FtD27KYXiQSiQwFEoIpBtk2AxsuVMCvAN+ovX8V8ElVvRD4ZHiPiDwcuBx4BPB04G21CqpvB14EXBi2Ml3IC4CDqvog4C3AH4S+tgOvAx4HPBZfsbQUxD8A3hKufzD0EYlEIkOFFRlo2wxsqFCJyNnAjwF/Vdv9LOCa8Poa4L/X9l+rql1VvQO4DXisiJwBbFHVz6mvAvmeJW3Kvj4IPCVYW08DrlfVA6p6ELgeeHo49kPh3KXXj0QikaGgDKZYD4tKRK4Wkb0i8rXavteLyHdE5Eth+9ETeDtrstEW1Z8Avw6Lsis+QFXvBQj/7gr7zwK+XTtvT9h3Vni9dP+iNqqa48t371ilrx3AoXDu0r4WISIvEpGbROSmffv2DXq/kUgkctyU4enrZFG9m74Xqs5bVPXisH10Pcd/tGyYUInIM4C9qvqFQZsss09X2X8sbVbra/FO1StV9RJVveS0044rEXIkEokcFSI+hdIg21qo6qeBAyd+1MfORlpUTwCeKSJ3AtcCPyQi7wO+G9x5hH/3hvP3AOfU2p8N3BP2n73M/kVtRCQBpvD/ISv1tR/YGs5d2lckEokMDUfh+ttZen/C9qIBL/FyEflKcA1uW/v0E8eGCZWqvlpVz1bV3fggiU+p6s8BHwbKKLwrgA+F1x8GLg+RfOfjgyb+I7gHZ0Tk8WGO6flL2pR9PSdcQ4FPAJeJyLbwH3AZ8Ilw7F/CuUuvH4lEIkPBUc5R7S+9P2G7coBLvB14IHAxcC/w5hN2MwMwjOuo3gRcJyIvAO4GfgpAVb8uItcB/wnkwMtUtcwL/FK8n7UNfCxsAFcB7xWR2/CW1OWhrwMi8kbgxnDeG1S1NH1/A7hWRH4H+GLoIxKJRIYGOcELflX1u9W1RP4S+MgJu9gADIVQqeoNwA3h9f3AU1Y473eB311m/03AI5fZ3yEI3TLHrgauXmb/7fiQ9UgkEhlaTuQaKRE5owxqA34C+Npq559ohkKoIpFIJDI4fsHv+iiViPw1cCl+LmsPfo3ppSJyMT6Y7E7gxetysWMkClUkEomMGEdVOHENVPVnltk9VFMeUagikUhkxIjZ0yORSCQy9GyW9EiDEIUqEolERgwRMFGoIpFIJDK8CHIK+f6iUEUikciIIQK2Ydc+cZMQhSoSiURGDSFaVJFIJBIZYkQwUagikUgkMsyI2egqTSePKFSRSCQyYogQLapIJBKJDDdxjioSiUQiQ4uIxKi/SCQSiQwxArJOuf5GgShUkUgkMnIIxsZgikgkEokMK3EdVSQSiUSGGYlCFYlEIpFhJ7r+IpFIJDK0iAg2jUIViUQikWFFQKJFFYlEIpFhJmamiEQikcjwIrEeVSQSiUSGGImuv0gkEokMNcIpFUxx6txpJBKJbBIkZKYYZFuzL5GrRWSviHyttm+7iFwvIreGf7ed0BtagyhUkUgkMmqEBb+DbAPwbuDpS/a9Cvikql4IfDK8P7ahetFba9u6Wh/R9ReJRCKjxjrOUanqp0Vk95LdzwIuDa+vAW4AfuMYL3FP2FZTTQucu9LBKFSRSCQycsiJrvD7AFW9F0BV7xWRXcfR1zdU9dGrnSAiX1zteBSqSCQSGTF8hd+BhWqniNxUe3+lql55Aoa1Et9/vOdEoYpEIpFRQwTTGPjxvV9VLznKK3xXRM4I1tQZwN6jbF+hqh0AEXkgsEdVuyJyKXAR8B5VPVSesxIxmCISiURGDu/6G2Q7Rj4MXBFeXwF8aB0G/XdAISIPAq4Czgc+MEjDaFFFIsfJvsPzAJy2ZWyDRxI5ZRAQuz6l6EXkr/GBEztFZA/wOuBNwHUi8gLgbuCn1uFSTlVzEfkJ4E9U9c/WmpsqiUIViawTM/MLAEyOtTd4JJHNjiDrGfX3Myscesq6XKBPJiI/g7fQfjzsSwdpGF1/kchxstSSOjAzz33Tcxs0msgpgYAxZqBtiPgFfNDE76rqHSJyPvC+QRpGiyoSWQdO2zLG3iBOLuzbd3g+ugMjJ4xRy/Wnqv8JvKL2/g68i3FNolBFIuvErqnx6nU5bxXFKnIiEBFMeuo8vk+dO41ENoh7D81ROAXg7O0TGzyayKZA1m+OahSIQhWJnABO2zLGbftmaFnBGp85xgG37p1BaolkwqFFuWUUcAqq/X0X7po80UOOjBIjVOZDRF4NfFxVB4rwW44oVJHICeJBp01yx/4ZsmBN5a5/TBXqz5lSnFbKIXrL3sOYIGdRtCJwVJkpNpo7gF8RkUcBXwY+BvyTqh4ctIMNu1MROUdE/kVEviEiXxeRXwn7V0wvLyKvFpHbROQWEXlabf9jROSr4dhbRfxvVhFpisjfhP2frydeFJErwjVuFZEravvPD+feGto2TsbnEdmczPQchXohypySOaVXeOFyunhT9WJWt6YcikNXvkDklETkhC/4XTdU9VpV/fmQ7+9PgQuAvxeRT4vIa0XksWv1sZF3kQO/qqoPAx4PvExEHs4K6eXDscuBR+BT0r9NRMoVb28HXgRcGLYyZf0LgIOq+iDgLcAfhL624xe1PQ54LPC6miD+AfCWcP2DoY9I5Ji46MwpHnTaJBfs7FtBC5ljIfdbr1Dma+/LfaWg9QolL6BwUKhSaBStCFUKpUG2YUJVv6iqv6+qTwaeAXwd+KW12m2YUKnqvap6c3g9A3wDOAufXv6acNo1wH8Pr58FXKuq3RDWeBvw2JCHaouqfk5VFXjPkjZlXx8EnhKsracB16vqgWB+Xg88PRz7oXDu0utHIgPzyg99jVd+6GuL9n3qjvv51zsPsHeuy/75HgfmMw51MuazgpluzsEF/3ohcyxkjvmsYK5XMNPLKVQryysSAUbGoloJVT2sqn+nqi9a69yhkNvgkns08HlWTi9/FvDvtWZ7wr4svF66v2zz7dBXLiLTwI76/iVtdgCHVDVfpq+lY34R3orj3HNXLKMSOcV44v/+F1yhPO5h/s/2R9/xWV7y3y7g6s/eySPPmuKRZ25hpleQOcWIkBohrU1MZYXitMCI4FRJw4NGlSri4pa9h3nIri0n+9YiQ4SIYNYphdIosOFCJSIT+GSF/7eqHhZZsbbWcgd0lf3H0ma1vhbv9GnyrwS45JJL4u/cCN/32/9Ea8xnhOnljh9/5On0QgTFoZkuU2Mp39o/hzXCRCvBBqEaSy2JNaTGRwimRjDiXxtxjKUWF5RqdObPIyeaUYn6Ww82VKhEJMWL1PtV9e/D7pXSy+8Bzqk1PxtfNXJPeL10f73NHhFJgCngQNh/6ZI2NwD7ga0ikgSrqt5XZBPQnZupXjfH1z96bv5wF+eUyx62i+luzg8+aAfTnZxnXnI2ew4skAQxmu3ktBuWRmIoFFLrSI0htVJZWmOpxVhvWVnjRcogiMCeA7MAVeg7wBlbx1caVmSzMULh6XVE5CJgNzXtqT37V2TDhCrMB12Fr/74x7VDZXr5N7E4vfyHgQ+IyB8DZ+KDJv5DVQsRmRGRx+Ndh88H/mxJX58DngN8SlVVRD4B/F4tgOIy4NXh2L+Ec69l/dLbR04iCx/5CxoXPALaW9BGGze2DbUNMBbkxH25b3zdZbzh+lv4njO2sHe2S6dwzHZyZjs59053mO1k1bnd3DHWsEy0UnZNNivRaje8O6dlDa6dkoaH0Viq5AiJVVIRHJCaFb0PkU3PCa/wu+6IyNX4GlRfp59pTIHhFSrgCcDzgK+KyJfCvtewQnp5Vf26iFwH/Cc+YvBlqlqEdi8F3g208TH6Hwv7rwLeKyK34S2py0NfB0TkjcCN4bw3qOqB8Po3gGtF5HeAL4Y+IkPMzHtej0kTXOanFpPtp+HmDmNcgXE5mrTA5mhjDPB/MpJ36PXmaGw7fV3H8qbX/gVv/dNf5Z7pDr3cMd8rmF7I2HNgnvm5Hhq+nmJgOrW0WxnWCO2epRmEqpEYFowwlxW0EksrMUx3clqJYSy1nDbeoGF9JJQ1Ui0g3nNgNma+OEUQc1SFE4eFx6vqw4+l4Ybdqar+G8vPCcEK6eVV9XeB311m/03AI5fZ32GFOiqqejVw9TL7b8eHrEdGgENXvgbbWrzUTecO46bvR/MMUxSYpIlrToAxYBvg8urc/J5bAEjOfMi6jekjX7mXBz9gktwp0wsZs52MXq+gt+Cv69QHUrhCMUZY6HnxLOezermjkRh6uWOWnEZiyArHllZK5pSJRgL4OS3nFCuCxrD1U45Rs6iAz4nIw0Ny2qNi5CQ5cmrR+/+u8y47oPH9z2b+7/6I+e/cy8LeQ4g1TJx1Gi7Lmbn7u/6cLWPADOnsHMlYCzO5jTTPsJNb0bSNNse9GxComzfZvrtxzQmS++8kO/1htNrHVlPqB372p5mZzzi0kNGwhtlOxj375ugt5MzP9tAQXy5GaDT9fe093KHd8BbVQlZUfZXzT4VT2qll61jKRCshKxxTrZQHTDRoWD+fFZ2ApxgiiBm5qL9r8GJ1H9DFGyqqqhet1TAKVWRk6Hz8yup1mTn6uzd9E2MNYg0mTegemiUdb6GFQwtHChQTW9E8QxotpD2OAJI2QZ0XrSTFtXxgRXb6wwA4NDvP1om1s54/4U2f4qmPPhOA337aQwF45DlbuW96gXYjYb5XYIyQBwGqC5VzSt4r6PYKCqcUDUs3WFWFU2Y7OUXuUFWsNUyOpWwda9BIDKeNuyBYCWOppZ36X9e37p2JKZZOFUZPqK4mTPfQn6MaiChUkaFGkgaa98D5v2tJvTXkspzuoRlm752hyAqmzt0K+Hkb18txWY4W4bvgCt9HeE2SonOHMVM7kKKHa45DkYNtMJdrlW/vrvtnQ1g4/Nvd09x90Jfu+L3XvZ3t5z+KhzzuQYvG+uyrP88FuyaYXsiYaKVYI8x3croLOSYx2MIBBnWKTQRjBFkSEFFmWe/lju5CRp65cH7fzbPvcJdGYphoWEzXC2AjhLcD3L5/ZlEmjMhmRLwre7S4W1U/fCwNo1BFhhpJUjAGNzeDJCnamcOmaXDxwZbCUWQ5RYioyzs9XJbT3OaDCqTRAmPRhTlksuEtqyRFxv2CWbUNNGmijTE0bdPrOmZ7jrmsYGc7oVfAH3/6dn7ggh2LxrVr9+kszHrx+9RX7sUYYdtUi+n5HtsnmrRTW7nx2hMNVJWsm+Ccok79ZLgRktSQ1KL9IIhUmNPKswLnlKR2/P45f92WNWQtv25ropFgBB/+LsKeA7NYIzFkfbMiIKO34PebIvIB4B/xrj9gyMPTI5FBUBei9KxFaSBjW2hsnSMZb9HcOsnWB58Tjhu0cLgsx/Vy0i1j2PFJtCgg73nrqdFCpnahadNnK0+8iHXHTyN3CoV6F5wqhYNP33WIA52Mre2UKz95KzMHFyhy5ZzHXMr+7xygPTnOxNYWzXbK2HiDHRPNatxjDUtihN27vGBaI0zP97h/yTxVufg3CdZQt5Zi3anilsmZVDitLK8SKxI2L1alobZ3em5RQcfIJkEEkpHLl93GC9RltX1DH54eiaxIccfNSNHr7zAGCi9aZnwLkvUwaQczuc1bTXkPtzBH0mx7956xSJJiSouq10HCF1tCyLq2JlHxC24BOoWykPucetZQleeYaCU0mgnNVsrhgwvkWUFrfAxrDUlqSVITFuR6dWin/pduN0TvlcJy+lTbh6I3Eu6f7dJMDDOdvFoEXN17uK61BtMUisKFaxnECO3Uh7CPpZbJhmUstVhD2LxIiayeyPPeQ3NAXCQ8qsgIrqNS1V841rZRqCLDi0m8CFXvLZKmEDJKmAk/L1XOW2nWgzxDXeFFyRhMu/8g1rSNa09VkX8d06RXKN1ugYjQK1zIWu7YP59RqA9omJ7PGGsl5OMpYoSs68PM02bCxNYWu7Y0aTcSJpoJ0wsZvbCYd7KVeDde7iqxajcSCqc0ExNcc23aqeXAnPeElOc5pzTbSWVRGSOMjzfYMdHgjK3+M3nuo/sJWe4+MFvNUZVzbHVL6r7pOVRhaYayew/NRbEaRYSRC6YQkWuAX1HVQ+H9NuDNqvqLa7WNQhUZOrL7voVJUj+ftGUnAEZdf/2TGLTR9vNLtoEaC65AXA7qfMxrCEEvALUp2AS1DVxjjG4opdHJiqqURlnIMHdKJ3eV66yRGCZaCbu2tJhsJdw/28MaoZEYmolhaqzBGVMtCqfM9wrO2zHG9okGhVMWQjTfQq+oFv+WllMjMVhjaCSGyVbi0yp188o1qKreWkv8OTsmmmwfb7BrS5Pd24+MRrQia1pRdUT6Na/iQuFRREZOqICLSpECUNWDIvLoQRpGoYoMJZq0fPqGsKm6Reue6kEQAKLuiHMA1CT+mE1Rk9DNHZ1CyQulm5eVd5XESKiw63PrZYXPXD4eVv/vGG+Qt1O2jzeZXujRSPxD4uxt7X7OPqdsbaeMpZbZbk7DGnqFo92wfjHvbA+be8vNGqFhTRAsqcSvRERohn7HGpbtwZraNtZgqunH9JV7prnozCnAl7Q3od1SSjeflSNjgsuprrvun+W8HVGsRon1cv2JyJ3ADP53Xa6ql6xLx0diRGRbWdk31AUcSIOiUEWGDm1Noi4HxsEkXpgq0RJc0iQPLrKFvF8tNzU+QEFVKYKbyyI0Ej9vUzhltufde7nzc1DWeHEqxWo+K+jkzpfYsMKUTZhqJYw3klB2QygUOrm3jqaaZdSdZedYCJ1XZb6ZVHNcRajsu3O8QeaUhaxgoVdU2SdKa8yKdwX28oL5XsGuLS3aqWWsYTlre5upZspUK2HXeJN26sPR79g/U7n8zljFKpIgUqX7T6uqwlrN0d2+3yfsjaHtI4CY9Q6meLKq7l/PDpfhzcBnReSD+CCKn2aZTEPLEYUqMnRo2oYiA5tSmBSnodqt89F4vU5WCc1C5rDGV8BNrXd/WfHCZMVnHE+MVC6xnlM6eb9SrhZCwwqZ85bVTC/HKXRyx0RNnKaaSVUjqk4zMTiFsdQw0UgQgV6ujKV9t4wJVs5sL2c+K5jPCqa7Oc5plbh2spUw08mZIGEhWFjt1M9zTY2lTDVTdo6lNJMywzokZnH29FU/U/WBE985OFdFHZbCtfSeIiPACIanq+p7ROQmfHFaAX5y0HRKUagiQ0dzfJKFTsdbTYXDKXTDz/5u4ZjLHL1c6Rbe+vHCJNiQfciIL5EholgnpAYcpXXjS7qr9oWsFC1ftBDmsyJYWxoKG/pce77AoZAaQ+b8uFqJoXD4yLswTzTRlP68l/STxqY2JbVCM5T2qELRwxxVu2GriL+yDEjpFmyGObFWYrx1pCsUSluFuw/M+gXR9OenSpEq+4rW1KhwVAt+dwaBKLky1NMrUeCfRESBdy45tq4EYYq5/iKbg3arxfTcAr1Cmcsc3cK7z3q5Mt3N6ObehVc4rayKpvXWhg3K4KvjKodDVou65eBrPvkveisxOLwbcCykImqG8hrWCFnIcJEV/lqpFSaDy69wYFOfGaKd9K0bg+BQDIutuomGoZP7RLIzvYL5zGB6Oa0ltYXO3TbGWGrInNK0JoxDSY13JapA5gYrTZ/VTqrckdoXKyCmXRo1ji7qb/8a805PUNV7QjX160Xkm6r66eMeY0BEblbV7z2ec6JQRYaWqfE2d+yfoee8628+K+iGcO/ygZs5rQSoKt/uw/fIxB/zZd/9Q91bRgYjMJZSWUe+sq6haX2WB29dKZlzFNZU10itMNlIKhFK/bQZaQiOSA3VnE/mtCrFYQUa1gdsFA1/fqdwdDLHgYWMzGkliK1gcT1gvEHm+uEPE42kstB6hWMh08qtCPDvd/lKNY8/b/uizzEvXX01ccoGUbjIELN+SWlV9Z7w714R+Qd89Yh1EyrgYSLylVWOC76o7YpEoYoMNakRcucj1qx4a6aphjSIDiHark5WOIyYyqUHfcvDiPjgBtWqxHsjWFYi0E5NlYqoUKWphqzwbsa6FZaGOlBlcII1hIW2QiIs6sOKnwdrWKnGmjcMbScsBEuqkzsy57zVZP11jAiTzWSR5eNQ1HlLbrbXL1dSVgYG+MQte2mFCMInPXAnhesHUBTq3ZvlZ2Ni3vXRZR2i/kRkHDCqOhNeXwa84bg7XsxDBzinWO1gFKrIUHP29glu3TuDCIw3LKqWiYZWIeQzvbyyfgqfBQkr/fkfE8QNqB7ehWORhVS650pBKY8VKqhCwyo2931ZEcbSUFpDjhQiK/48E2Z9THdxliUAACAASURBVOewD5nXBuRQ/gbeobnPYTg5yVhqyENkYO60EtXUCInxEwid3M+XZU5ZyB3zWcF353rkwQprpzZYhbLIyvrCtw9VmdXzYH0WtRj1dLTm4yMlYqpMK8fJA4B/CMsaEuADqvrx9ei4RFXvOt4+olBFhp4Ld01yRwid9giKd181k77FkxW6yIoq56pKsWrV1imVD2u/nqkvVEl40C8N4W5YX/MpDYENVvyDfyKxlXVlevOLxi1FD+kt+De9BT+nIMYvTs47oA4zf5DTGuNoewpnUxYyFwRXF43VqV+EXGqMNQQXqMM57xYtaSa+ttVkw5Ja8a5NkcpiK4XMGnBqFglbZEQQ1sWiCoViH3XcHZ1golBFRoLzd0760OraA9xnoCjTAlmMuCNCrat5K9N/IDvVKqS9dN0J3hJLTT+pa+b6V7JOq7moMvHrWCI4BJt5gZKi58uFANgEyXtenFzhFx3jX6POHwPIM4w6tMgQmzJhE1xrChDmcz+3Nud0kduwjqsFR2TO0csdbadkzlTWVjknB1SuRSsCTgBXuTMjo4MgIxeefjysKFQismqURiBT1a+u43gikRU5a5vPSXf/zLzPIuEUcCFruJK6vhBBX6RKC8kgiyyVVmKqCrlWCOHfQio+cW3bgrMpoopkHSRMDkl3zrvzxGCsL0OCc5je3OLMGGKq81CHZAs+zVOeIUWvX3Kk6KFmwS9sTpu+f2MZb2+lsC3A0SmgG9ya1oBxpUvT0s0LDneyavHwgdkejVA6ZD4rKgsrDaHu3YIq7L4ZAtZvvPsgzcQsyrz+0AdsOTn/sZGjZzRz/b0ceH+ZmeJoWM2i+lfgRlh1tvV8YPfRXjQSOR52TI5xYGYeawXwIdwSPF9LE69aKYVIgqvLzzuVUXvNRNjSMJXLz7oMybo+r6BJsN0ZsA1wuRerIvNCYi1qwtenwAtR0fMZ3kWQIvPZNJLUF310hT/uCjTPcFkPSRtorwhWl7fKTHMspIdKEZNgjaVhU5w6CkM1NzfdyX06qOD+q+cVBO/S7OWOTuHCmq2CsdRiQmXhiWYCeGGbbPoHnhd4CVk+ojtwuBnJXH+nAzeKyM34ar+fUNUj3QTLsJpQ3aiqP7RaYxH51OBjjETWj/I5WuboMxJqPAkYlWqBr8HPOzWCu0trgRJlyHjLSiUiZWLbyoVnEi9SeQ/Ju37eyeX4RzpQ9PqWUyhDIkXm25gErYmU5hnaXfDCVXe3ucJnfgefIV4MqPXuQldggiCWllShSrfwARW9wlUuwcJplS+wtKjKGlUAncKRql9r5iMjrc8RGBZMm1pGjzhtNeSI+KKiI4Sq/paI/C98ZOEvAH8uItcBV6nqt1Zru6JQrSVSg54TiZwItk74ZLQHZuZJDagKJCHwAPXv8Q/d8dTUxKrfRxmdZzrTlRDhch9RVfR8BvaQZ9DO7qvqXIFPgoszfdFxzgdIlISkuJJ3QZ0XqYU5NK/V2Mp9VWJ1he+njhgoMiTMfxlpVOmOUuNrUU00kpBFwzE1li4qqFhmtCi8+di3KoMwNUMQyFhqSY0Xr/Iziho1IsjozS2qqorIfcB9QA5sAz4oIter6q+v1G7NYAoReQLwJVWdE5GfA74X+NP1CDmMRI6X7ZNesO6b9hnCMxTFP2xNiOJrhtDxBOetoiBGGgRJunNI3kGybmVZYSzCXHUd7XW81eMc0gw1soxFkrxvUQEaxIdaIIUWXog071WCpK5AyRBjg1AVaJZhjEHSBlL0MAvTqMtxQNJqkpggOOHfqVZCM/FBE53CR/+V6Z7KBb0miFNqDe3U0gpzVVOtpBIov1DZW1H1vIiRYUZGTqhE5BXAFcB+4K+AX1PVTEQMcCtw7EIFvB14lIg8KnR0FfAe4EnHO/BIZL3pR/n5dU5l1VtrBCmcFymXe8ExBilyTDZfzUtpr+OFqk5Yr+LmZ3w2AFdAkiLWehegMd7FB1Vb7Xa8AMEicaooz4PKshJrq/41z/zcWIgQTEJEYh4yYEy1Eopwn2CZKeemBOYzIVWt1nh5C6oMV/eiVIpUWRHYmjA/H4YnIrGo4pCjIyZUwE58ItpFRo6qOhF5xmoNBxGqPJhrz8JbUleJyBXHMdhIZN05fWqcew/NhQW13pICwpwLiCu8K60334+4yP3ckSwcxi1460l7wX0XXHpirHfZdeb9/FKzjZs7XFlV0mj5hZel8LjCW18Lc168SsvKWD+ncMQEeOYFrwhWlSsq15sUPciCC9HlNJMUa8rUUaaqd5UVylTLVevIJhpapZYqo/tKK6rMi2hrk1DNkNWidIuK9CP/IkOKMHIWFXD+UpESkfeq6vNU9RurNRxEqGZE5NXAzwH/TUQsMFqzeJFTgqXiZMO8jEG9Wy/vgKoXgDJk3OWoc9WalCoEKQQ7uM4cmmVI6v/kNfPCo1lWXVeCEGmvg+Y9NMtw8zPV+WL9cepCVc51GYO6IGTNIHrlOeUYewvY6XswaRvbGMM2J3xBxsRUeQ8XakUgZ3s5zXA/9dIn9aS5ZeL2Zm3f2dsnqiKL0ZIadoQRjHh5RP1N0JLHDNJwEKH6H8BzgReo6n0ici7wv496iJHICSat/cCsZ5iQzM9LSYjk01CMUVzu56mgLxyNVmVVeTegFzHvyvNl7iVtoEXRFyDwwlZG8pVBF7kPQSf0X7avKNuG8zXL/DxWmKfSrIskQVEKX6GYrEPicmzaxiYNwGCd0kp8Fo3ZrACSaj1UaTmV81pSRT6GoItwvFyjFgVqNFBA7WjkawiGzmuAtogcLncDPWCgkiKr3mlQvPep6lPLfap6N36OKhIZKqzpR6xZI4grkCJHsnkvUq4I8eup/5Ln4fzmeD+IwhWQNPrReSGirxKvcNyLlPHuvDBfhbFozVKi0Q+6qCyu0rVo+laWtLw4VJkGnKsCN7RMmGsb2OwgrtFGm74kRyNRGtaiiUGNr2U1ngqdwufzM7Uf3WWZeiugCM2QGgpg11QUp5FDRieYQlV/H/h9Efl9VX31sfSxqlCpaiEi8yIyparTxzTKSOQk4ddTLbNCXQxqE59ySQzYBEyCmgRVh5TRf3nPu/GSDOaKSmjKDBJlwIT2Ol6k2uP9KD9j/evSOgrnShAa7XX8HFee+Xks5zBJCs1WELtG/1plwAb4NVjNFnTnIG36tVVF5kU474WEtwrGYm1CmrRoh/U11YJkQEPWeKWfvxCkCvOPjCAjIlQi8lBV/Sbwt8tlPFLVm9fqYxDbsQN8VUSuh368rqq+4mgGG4mcaKbG2wDMzi+EeAnTT2VEcJWIqdZH+Z0OsT0vYGXCWIBmC7rB+kmA9jgSAiy021nstiP14lOGn9ej+4IAifHi5PKMotNDC4dJc5K6S7AM3nBFEMyQ1aK8ljrvigzh9ZQRjP4C4Br+/LTlzzVJPxGuGIzYkFIKXFwtNeKMjkUFvBJ4EfDmZY4pvjT9qgwiVP8nbJHISKCE7OMiPquDWebPvHSdqKCNMf/Az4MgZBZsA0n87zJxBdpo9Rfo1lyDmmVVYIVYi2Yhh18IWdf5w9UlNfcClc0t0Ln/MLbVoL0jp73ttL7bb1HGirCYuCaKflFy4efY8k51LZImalM0baFFVgmUtkLl3tpDrUzAGxltRiU8XVVfFP598rH2saZQqeo1ItIGzlXVW471QpHIycapIkYQuyRItUwUW74V49MdJf6YmsQHX0BIfQTSbCNtP5ejeValPHKzhyqh0iIsFDYFLrgBy6ALLYKVhLeAeofn6NxxH8mjH0Q7aYT5rsbiqq010ZIgWPUACzd32LsSAWmP+1D5IsO1JlFpeqEyCS5YUrWYxljaYzMwIkJVIiIvwyelPRTebwN+RlXftlbbNe9URH4c+BLw8fD+YhH58PENORJZf2bnF5idX/ACQ/+xrCELRfXaNvz8VM0tWBH2qUkgSdGk5YMtGm1ImpA0vSA0236BbqPlw8pTH1Ah1np3X7l2Kq2FmweSVpPm1kmmHngmza3B4jF2sfVUBmU022EOq9lfq5Vni+bGvID1QpRiXpu36lcuLnP5TYy1mRjzLtKxdmvdPvvISUakH5Cz1jY8vLAUKYCQRf2FgzQcxPX3euCxwA2h8y+JyPlHP8ZI5MRSPoBn6mJVzlVRCpb0S7uLzwRhQqZzMQma9tMhOZv260eVKZJCWLv0Fnzhw6ThgyToLxaWWpQg1AIq8gzJMxKgMRUEKgQ+iLHeomqP91M4lS694LoU68B1+8eb7X5C2zJS0Fg/NvGiJ3kXsemiwAqIIrUZGBXXXw0jIlJmTA9R5QOVKR40M8W0LHYVDJSaPRLZCCbH2ix0+kKhtcCBelEBEZAQMYf6iDhMUsV0axIe5mUJj5LcJ6yVoucFMYS3a7PtM1GUJcKN6UfxldZNCJbQznyVjUJSX9NKmi2k0a4sPZekVXSiqPNzaCZZNBbNw7quMvAj73kLLPNzVAZITzt3nT7ZyPAg61Lh9yTzCeA6EXkHXkNeQvDUrcUgQvU1EXkuYEXkQuAVwGePdaSRyMliuUo3ZSFFCCKl7siTyl+qpb6poLUffpL3QARNvQUnRYgatI2q9IKZ2oGGtEymPe7nqIIb0BiDSxq1NVUGabS8Cy/UsCr70xC5p67w/UN/zC73QRuhD99XGSGoISpw+a94XcgB2q1oYY0Uo5lC6TeAFwMvxd/BP+GT067JIEL1fwG/CXSBD+BV8Y3HNMxI5CRRlsRYtI/w/Q6WlE9QGwoahnBuNbYKtlAx3pIR42dzg7BVUYLgXW8m8QEO6lBte0tmYdpbR70Fn/zWWDRtozbx40jbmKJXZWU37XGwKZqkPnweUJtW4/BZY2XxuIqkSgMladOPp5xfs3bFB1lnYSGsKetbmrPzCxiR6BIcGUYqPB2oks9eBfwb/ut4i6oWazQDBhOqH1PV38SLFQAi8lPA3x7LYCORk8FYu8Xs/EL1vpQsa6QvUkUIK89C1omE2rqrvlhBMK4EwGdMx3mRoGm8u09dP/OFKsX4DkQdkjZxZYBGUhMBk0Dhq/hKw1cTPmLOQbVv1dnEZ6kQ489VVxPHsLaKxYt8kzMfsuLnU87VlWIeowBHj6XzjsOOiFwKXAPcif/LPkdErlDVT6/VdpA7fTVHitJy+yKRoaRuV8kSK0uWc/2pA6klhuXIiWsBcARrCjRUBRbtpz1CXRV84RcZ18Sgsnz80lu1fp5KTbI4fL7+upZhAxVUTSWMVbRfYCWR6iyEwI8lFudS6zMy5IxQCqUabwYuK5c5iciDgb9mgMS0KwqViPwI8KPAWSLy1tqhLZSLQTYpIvJ04E8BC/yVqr5pg4cUOQaWRgEegbH+YZ80jtjvhcRWFXPL1ENVqXabgk19xKDpu9m0dBuapHKtmbwLgPRCYpfyAZM0/TxXWXhxqVgYuzibRomUFpzz6aDC+BrbTl/18+jMz1XXdqqL7qlkdn6h+twiQ846WcEn8XmX1tfiqup/ichAlThWs6juAW4Cngl8obZ/BvifxzLKUSCETP4F8MPAHuBGEfmwqv7nxo4scqxMhgfv/EKH0r4qF/Ui4EKSV1FXLfhFDHmttHtJoYABCbvVpIhA4ZTcKWDL1BiAUig4DQEW6daq1HvaqCWCVVeVrK8EyARXYSma9cwSLoci9wEUdctrFTrzc0fsK3P+OWRRjsQoVqPA+lhUJ/l5d1OYo3pveP+zLNaWFVlRqFT1y8CXReQD4bxTJTPFY4HbVPV2ABG5FngWEIVqxCkDBSr3V/DxO3w5ELR080klUGVZjCJMFxnxD/gyBZFTJS/8uYV6waov5VDVqg+H0Ml9rj2bpr6YI2G9V7lWyuVgywCKdLFAlWJURf05EFebV/P0Dt6HNsaPdFcGV6QGd2GZ82/ZRL6RoWed1lGdzOfdS4GX4SPHBfg0sGZWChhsjurpwB/hF2adLyIXA29Q1Wce21iHnrOAb9fe7wEet/QkEXkRPtEi554b16mMEq22X2fVd+upN2LKB3ZwwS0xpipDKZFgMOGFqVCfW7AodWSJdeMzQ4Cg+DANYT5zWPGViBPBL8oVg+JTHzmEzPnzBX9N1PkAkLAAucq2EazD7pwv1ohdZg1lPTBEFWu8GC/N+Ve6Aw/NzsfM6sPO4EK1U0Ruqr2/UlXLOlADPe/WA1XtAn8ctqPiWDNT7D7aC40Qy/24PGKmOfxHXwlwySWXxJnoEaPd6kcFrpSkVWuCVS8hoqrebQZkri9QhWoVpFDUrK5yDigxgpaCtWR+obSocu0HOvSCGqYGktoApfCBE34soUpx2l8ojElwSLVmrIy3CCuuyPXI3GmF9q2/yPCjyNFkwN+vqpescGyg593xICJfXa1PVb1orT6ONTPFZmYPcE7t/dn4+brIJmNirB3mrRZHvdW/UUutKvDuQEMQq1oBwnrghSqUC0Q05MbwuiPeqgpCVigkZdl7BMXPi7klwuEQFmVtUwVc3/1XCheEtWDl+LV6baTv1tTaXdbn4QqnlbBOz3khL8unRIYJXa9IzZPxvHvG8XYQM1McyY3AhSGf4XeAy4HnbuyQIieKpeut8tKVFx7eDm991KsH50FIvLuvDJhQcgeOUrQWXyc1QUaM1xgrkDkhNYrD4FRJRFkolMxpFQBoUZx6YSwX9+IKb0WF1EparuEyfh6qNTbOQqeDCe5D8KKX0xdk/9pfo7TcABpWKNRbcSX1zycGWQwP62T2nPDnnareVb4WkfOAC1X1n0NVjoEWgx1tZoq/ZpNnplDVXERejr9PC1ytql/f4GFFTiCzmatcYY7S9dYv5S4iiyyr0sVXuCAqeGHqiwJ0cq9UaXDZhTScKIoVoVco1niLKgvVgFPjr5sXWglkGb2Rq5CGsHcR6adWsomfkxJDc2KqGuNyKZH2HZ6vXpcRjaXIli7CXuEvaUSwRFfgsFLOJR53PyfxeSciL8TP628HHoi33t4BPGWttoPUo5rHC9VvrnXuZkFVPwp8dKPHETk5lBbFEVPTZXJa7VtX/r0XqfJ9UbekFs1TBWtIBGc0ZHEXMCGZk/NBDbn2Ba5Qf53CKUW4phNv5VibYmr1skqR0nId1gCYEEpfDwIRgljRD//rFcoZW8c5NDuPisSowCFE18f1dzKfdy/Dxzt8Plz3VhHZNUjDNYVKRC4BXgPsrp8/yARYJDIK7JvLscZbEc1EsCKkRmiId/l5N53/uVm6++qWVKF9V6AqdHNHJ3eVWFkRUusFSzV8iYL1lDlvzZSCYY23tsooQiO+varSTAwN26DVmvRJaqFaZ5VJQnOAexWC+OHFySnsbFsOdAq/CFj7EY37Ds97q1KXDzaJbBzrZVGdZLqq2iutdBFJGNCDOYjr7/3ArwFfpf+jMhJZd/YcmK3+wLq54tBQTwoMwoW7/MLcO/bPYMOM/7nbJ477uo85Zyufum0frcQACQ1LtcaoDGBIjRcrDRZPKVIuBFWAP5YVZaSgVgIGPrKimzuaiUGcD6gow9Z7hVaRebnrz3Hlzu93AtYBuXcINtMGQm9RyHmCozs3Q3N8csX7NFJG9/X3lQJ02ljCvnmfcEZEsCtEPsfFwENC7QfFCPGvIvIaoC0iPwz8MvCPgzQcRKj2qWqs6Bs5aeQr/FS8ff/MES6ouw/MrotY3X5gnolmwq7xBlPNFJca2rVvR2L8w7uT+/mnck6qtKKWTuUUqmTF4sisQpXMGVJjvGAFAa7S9YUfl62gEjZYWrkqnRychcwVgKVhG6RGqpDzMiltd3YaXEFzy/Yj7rFpBcXfQ9m/qjLddWxvW3aNJRzoFBj1rkld5saSaFkNDevl+juJvAp4Ad7oeTHe3bhuZT5eJyJ/BXwSH1ABgKr+/dGPMxJZnr3TcyGzORRLIucMIflCOLd0oa2XO+pH3/FZHnXOVqbG0mp9kzUpYENZEKncc9CP4HPgzZ0wDgWwkBUE68zRzbUKsiB3OIXUKJlzpMaQWsEaqcSulyuqblFW89IVWV6jUGUhV6am+mmYurPTlGVIVmIiFJTcIl202URFmO4WUCjTnYKplmV7yzLddcHik2qNWRlsgWq1sHg16y1yYgmLE0aNZwHvUdW/PNqGgwjVLwAPBVJq88lAFKpNSJleqNU+ee6d+2d8NJq3MPw8TbHCr8XqgbnO3DvdoZs7zt3aJiv0iPByF+amTJi3cvhULWrKnHleQLyfzoee2zIprQhZUUYB+vPSoG6F6y8ILm85q/l00kRoWL9Vn89yH4BbXNYn23f3spV9260WtFp0Zw5B0mSbzTmgCZMNb8VNd12IdpQQ/efbLarhFRkKRs+g4pnAn4jIp4FrgU+o6kB/UIMI1aNU9XuOZ3SR0aAUqZPFgZl+uHRSsxhKSneYwQc0lFYFwNnr4O4r+ehLfoAnvfkGOGML81nBROHInAsh5oaktq7Iil9rVI63XFOlqmQOsJC7YJWJeDdhSFmUGlmUqTxzjqZdtIyXdriQNdAIbRIrtKxUOQOBI9IbSdGralpJKO6Y33sryRkXLnvPzcmt5RU5A+jOzaAmYUfLMJMpqemvHUsEv3arTJoLi0LhIxvDqAVTqOovhGzpP4Jfq/U2EbleVX9prbaDCNW/i8jDY/bwzc/JtKJKzDLGQSkEqdEjvoxnbRs/ssE68K+/eikA7/z8XWSF4zuHOzz+vO3cvn+GpjW0k350QUKwkpzSsKV4CJ1giVmjgJ+Dsl2pIvugvz7LiJAaQ+Yc1hhaiVlkKSUGJhrWp08ywth9Xwd19M5c/jdjY9vp9PbvwahDsoWjzqzdHJ/0yWxtgy1AYccWz3/VypD0RS6yUVQW/IihqpmIfAz/m7SNdweui1D9IHCFiNyBn6MK6c5ieHrk2Dk0O19lIl8qVoLf17CyKEPEyeDFjztv0fsLdvp5mOm5hWquxqCo9OduCHn10hDhZ1TohbByI4vn0ozg56ZMX0gK149ZsGEurmFN5fIzReYLIwKNe8NazAuOTN22yC2nbtUKvysh6iiaEz6gIiS/rURqwLVakZPDqOlUqHt1OfBkfO7YvwJ+epC2g2ZPj0TWnSJkgFiaELVM/moExPhINcvinHQnm3q+u+7stB8fVJaLMQkES9Dn6jOoOsZSW1lO4EXKGqrwb6dCOzUYhMQQLCsv0i3XhRAgbx70+DXHmO7afVz3WFplSW8e19qySJiWiyKMbBx+HdWIKRX8PH5u6sUhk/rADJKZ4q61zolEjhZXiZR3l3n3mVYLbEtK8SrnZu6b9gUAT586MS7AgSgDF8RH9tXdbD7xq/83MTCWWjKjdIt+FGNWKB3pR/apWr/gV33ARDMxNMUhvR6NyYEW7q8bjZ1nA1B8+6v+34nTwCT0DtxDY/uZJ3UskdUZNZlS1ctXOy4in1PV71/u2Gql6G9W1e9do+M1z4lE6hya9QEUZZh1mX28zLF3hC9QtXKLrUuZuPVEQ0Y+6b8uLcR+4cUjm6VWmM+82FkRHF6kfCmOYEmq82XqNwh7TpgL+/ZXKcZ3bNg4IiszasEUA3BkgsrAahbVw0TkK6scFyCG/kSOiaXP72ouKqQVKjM/lDSst7SsOTHh6UdDc8t2uocP+DfqKqPKGEvu/FopR1+wSsupdP9lDuaz0hXoBar+gSxkjoOFrMtC5uPFnvM9i8uLRIaG0fP8rcmKd7SaUD10gI6LtU+JRBZTFh00+D+gcl1SmSi1n+uuv55pvOGDC+yQ5J1rbtlOb3p/PyGs8/NU5Rqwwim5Uzq5Iyv8At/5zFWZykvxaiWWrNAq3P1Ax3+lHvqALRt2b5HhR0OKrlOFFYUqzk1F1pty3VQ9oWaZPqjcX+a7K9MPlWuNnAYLyywxPzaA3sH7qnkpgSrfHtTn3srX/h7K+/FZ0ZXU+pD08mGTO+X++YzJpuXis2L4d2RtNqHrb8Uv9kBFqyKR9WLpd8tKqGqL/yvNnDLb82HWhevXcyopHFywa2NdYo1tp1evy3RCKoaszJxeKN3cW1PzmV843M2LRbWfrFH2z/c4Y6LJbC/nSQ/cuVG3c9wUd30ZAHveozZ4JKcO5Y+6UWK59bgicqmq3hDePm+ltlGoIiecegaKMotDNf2k/Wq29WzQfuEstXx4J224R0WZ7+6+6TnyoiBzPg/f/vke81nBwYWsKo7onNJMLKD89EWbJ4IuCtTG4EYu7o/rROS9wB/iAyf+ELgE+H4AVf3aSg3X/PqLyMtFZNs6DTRyirDv8Hy1lZRfqzLCr9xgcfkJU6vfZI23ukzwCpSlPoaJW/fO0CuUTuFFaiFzzGcFM72C+cxv3VDxt5sXPO97z97gEUc2A/Xv0GrbEPE44Bzgs8CNwD3AEwZpOIhFdTpwo4jcDFyNTyQ4XLcfGWqKmhVVlnmvV83VWgDFQuYYb9jqC1YGGWROefjpwxdgcMvewxiE2Z5jrlcw08vp5o798xnT3Yx9h7v8+7fu5/EP9CHev/20QWKUIpHVGdEFvxmwgE+d1ALuUB0s3ckgC35/S0T+F3AZPpP6n4vIdcBVqvqtYx9zZDNTlncvBalMqFov517iAw68UDUS4XvOGI1VD7funaFw0A3BEgcWMrqFY7abM93NWOgV3PC1+4AoUJH1pV6kc4S4EfgQ8H3ADuCdIvIcVX3OWg0HmqNSVRWR+4D7gBzYBnwwZL799WMfd2SzUreYfCFE/97WQrNLukU/1HYUROrWvTPVa1/F10f2OQ3h6E7p5Q5rpEp2G4msLyMZnv4CVb0pvL4PeJaIVAEUIrJNVQ8u13BNoRKRVwBXAPvxSQR/LWTANcCtQBSqSMUtew+ThvIQZRh6+YUyCEX4gpWZwodxzmkQHH6NVy/3ItUNApUXjqxwfOg/vg3AK5/4wA0eaWQzMoquv5pI1fe9t/b2k8CymY4Gsah2Aj+5dF2VqjoRecbRDDRyalDPKOEtjvCa/kJeYGTXC5Uilbv+Gqmyiq8JSprMJgAAFPhJREFUc2pihBv+55M2cpiRzYxyRHHPTcCxr6NS1deucuwbxzqiyOakcF6QynR9ZfXbMpS2l/t/H3POaInUbftm+hGKriyY2Hf7leJrRDBGeP4Tz9/A0UY2OyfLohKR1wMvBPaFXa9R1Y+eoMsdUwqlSOSoKbMrmFpCvjIrA8D3nTuaKx3KZ0KncJXYdosiVAH2GBFaieHemViuPXJiUais+JPAW1T1j07WxZYjClXkuHnvzXvYOZaSWkPTGqY7+aJ8dqOcdaFkIe+X5egWBYULi3hrmdLHUsN3DnfZ3m7EtVKRE0vt724TEVMoRU4OBzsZTWvoFo5nPvz0tRuMEIULofTF8hFXVoRzp1qbQpgjw42iR+P62yki9UCGK1X1yqO43MtF5PnATcCvrhSZNwgi8oPAhar6LhE5DZhQ1TvC4aes1C4KVeS4Ka2H9968Z9H7zYSIT+nUzXx0X/lr1hrB1hYlP/XC0zZymJFTiKNYRrVfVS9Z6aCI/DM+scNSfhN4O/BGvLfxjcCbgV88qoH2r/M6fMqkhwDvAlLgfYTsFKp6YKW2Uagi68ZmFKgSDVFWRoTUCEaoQuxhc7g3I6PDegZTqOpTBzlPRP4S+MhxXOongEcDN4fr3iMiA61PiUIViazCngOzR+wzIlWwyLAmy41sck7SHJWInKGq94a3PwGsmDh2AHoheYSGvscHbRiFKhJZhTIXYTsxtBMAy775XnV8opFw0ZnDn00jsrk4iVF/fygiF4dL3gm8+Dj6uk5E3glsFZEX4l2IfzlIwyhUkcgaOFWs8Q8GVW9RlW6/br75Vl1Ghp+TtY5KVVesEXUMff2RiPwwcBg/T/VaVb1+kLZRqCKRAbEiGANjqa0Sgo7awuXIJkEVN2Lh6SJyPvCZUpxEpC0iu1X1zrXaRqGKRFZgz4FZDD4tUmk3GWB7ywJw9vaNrTQcOXVRjirqb1j4W+AHau+LsO/71moYhSoSGYB6hN9Z2waeA45EThijlpQWSFS1muBV1Z6INAZqeOLGFImMJnfs92U8rJGqBHaZu/CMrVGkIhuPr0c1cvOj+0Tkmar6YQAReRa+KseaRKGKRPDiVGhfkAQf/psmhtKYOn0qilRkOBhR199LgPeLyJ/jv2LfBp4/SMMoVJFIjXJ+OlpQkWFn1Fx/oSL840VkAhBVnVmrTUkUqsgpT+Xqk/6v1BELqIqcYugIVfgVkZ9T1feJyCuX7AdAVf94rT6iUEUiAQW6uVauvoefvmVDxxOJrMhoZU8v3RLHXM77/2/v3oP0qMo8jn9/uZDghUsuYkjAhAV3iyBCGAPUeskCQogueEOzu7UiWJUScdVyWTWiUovAbsRavK1CVigFoSKoLHiNQUGrXLlErgkQGUDXICoYQFcUmcyzf/Tzmp7hnZk3M+/M2/3O71PVNT2nzzl9zvsm9Ux3nz6nI4FK0vnA3wJ/Au4HTomIx/PYauCtFEMX3xkR6zP9MODzwK7AN4F35XQcM4BLgcOA3wBvaozLl3Qy8ME87TkR8YVMXwSsA2ZRzDv1jzkCRcAngBXAk8BbIuLWcfworAL6+nPS2Zx49q/2coCyagvqE6gi4qL8+a/D5ZO0OiL+rdmxTs1UtgE4KCIOBn4CrAaQdCCwElgMLAc+I2lqlvkssAo4ILflmf5W4LGI2B+4AFiTdc0CzgIOB5YCZ0lqrNq3hmIxsAOAx7IOgONL9a/Kc1oXu/dXv2XqlJx0NmDalCGXxDGrjAj4U19/S1uNnDTUgY4Eqoj4TkQ0lkG9EWhMu30isC4inso1SnqBpZLmAbtFxI8iIiiuoF5TKvOF3P8ycHReGR0HbIiIbbl+ygZgeR47KvOSZct1XRqFGynmpJrX/k/AqqQxynf/uc9l/7mjvjthNmGCYqmZVrYaGfKvxCrM/Xwq8K3cn08xZLFha6bNz/3B6QPKZPB7Apg9TF2zgcdLgbJpXU2ODSBplaSNkjY+8sgjLXXUqqdYqTc44HkOUFYj+YyqywLVkI0dt2dUwy3GFRHXZJ4zgT7g8kaxJvljmPTRlBlNXc9MLFbIXAvQ09NTq38Nk9n//PQ3AEyfMoXpU3PRQ397VjN1eka1EyZ+KfqRFuPKgQ6vBo7O23lQXMHsU8q2APhFpi9okl4us1XSNGB3YFumLxtU5gaKN6H3kDQtr6qa1dXsPFZTX7/nV0wVTJ86hWdN33ETYXs/HDLfS3RY/US9Rv216qqhDnRq1N9y4H3AKyLiydKha4ErJP0HsDfFoIabI2K7pN9JOgK4ieJt5k+VypwM/Ah4A/C9HA24HjivNIDiWGB1Hrs+867LsteU6nqHpHUUgzCeKC0aZjVy2a1b6Y9gisSeu04fcGz6lCJYeeZzq7O6BSpJMykGri0GZjbSI+LU/HneUGU79R7Vp4EZwIZ86evGiHhbRGyWdCVwN8UtwdMjYnuWOY0dw9O/xY7nWhcDl0nqpbiSWgkQEdskfQS4JfOdHRHbcv99wDpJ5wC3ZR1QDHtfQTGI40nglHZ33CZe40Xe6cCTT/dzzAGzO90kszHpj6jjWmiXAfdSDHQ7G/gH4J5WCipq8nZzlfX09MTGjRs73QwzqwFJP46InrHUsdf+i2Pl+etayvvJ1x085vO1g6TbIuJQSXdGxMGSpgPrI+Kokcp6Zgozs5qp6TOqp/Pn45IOAn4JLGyloAOVmVkN1WWuv5K1OWbgQxTjAZ4DfLiVgg5UZmY103jht04i4nO5+31gv50p60BlZlYzjSmU6kTSXsB5wN4RcXxOmXdkRFw8QtFKzExhZmY7oXjht7+lrUI+D6ynePUIinle391KQQcqM7O6iYmZ60/SSZI2S+qX1DPo2GpJvZK2SDquhermRMSVQH/RheijWCVjRL71Z2ZWMxM4hdIm4HXAReXEQStd7A1cJ+mFpfdem/m9pNnkpGU5gcMTrTTCgcrMrGYioG8CAlVE3AM7VuMt+fNKF8CDOeHCUooZgobyHorRfn8h6YfAXIoZgkbkQGVmVjMVmJR2PsUSTQ1DrjTREBG3SnoF8JcUE9BuiYinhyvT4EBlZlYzEbEzo/7mSCpPnbM2V38AWlvpoomWV5ooneck4Ns5Vd4HgSWSzmllFXUHKjOzGtqJK6pHh5tCaaSVLoYwmpUmPhQRV0l6KcV8fx+jWEX98JFO5lF/ZmY105hCqYMLJ14LrJQ0Q9IicqWLEco0Blq8CvhsXq3t0srJHKjMzGoo+qOlbSwkvVbSVuBI4Bu5fBIRsRlorHTxbQaudDGUhyRdBLwR+KakGbQYg3zrz8ysZiKgf2JG/V0NXD3EsXOBc3eiujcCy4GPRcTjkuYB/9JKQQcqM7PaCeq2RFMukvvV0u8PAy0tTOtAZWZWNwHbazbX31g4UJmZ1UwAMXnilAOVmVkd1e3W31g4UJmZ1c0EDaaoCgcqM7PaGfvQ8zpxoDIzq5kI2L598jykcqAyM6shX1GZmVmlOVCZmVllRYQHU5iZWbV5eLqZmVWaX/g1M7PKCk+hZGZmlRYeTGFmZpUW9PsZlZmZVVUxKa0DlZmZVZVv/ZmZWdX5PSozM6usiKDfc/2ZmVmV+YrKzMwqLfq3d7oJE8aBysysbiIcqMzMrLqCyRWopnS6AWZmtpMi6H/6Ty1tYyHpJEmbJfVL6imlL5T0B0m353bhmPs0DF9RmZnVzcTd+tsEvA64qMmx+yPikIloREevqCSdISkkzSmlrZbUK2mLpONK6YdJuiuPfVKSMn2GpC9l+k2SFpbKnCzpvtxOLqUvyrz3ZdldMl1Zd6+kOyUtmYjPwcxsZ0X/9pa2MZ0j4p6I2NKmJo9axwKVpH2AVwL/W0o7EFgJLAaWA5+RNDUPfxZYBRyQ2/JMfyvwWETsD1wArMm6ZgFnAYcDS4GzJO2ZZdYAF0TEAcBjWQfA8aX6V+U5zcwqpfGMqsVANUfSxtK2qk3NWCTpNknfl/SyNtXZVCevqC4A3ksxbVXDicC6iHgqIh4EeoGlkuYBu0XEj6JYLexS4DWlMl/I/S8DR+fV1nHAhojYFhGPARuA5XnsqMxLli3XdWkUbgT2yHObmVVH7NQV1aMR0VPa1parknSdpE1NthOHacHDwL4RcSjwHuAKSbuNV3c78oxK0gnAQxFxR97Ba5gP3Fj6fWumPZ37g9MbZX4OEBF9kp4AZpfTB5WZDTweEX3D1TXo2MNN+rCK4qqLfffdd8Q+m5m1T9DfpmdUEXHMKMo8BTyV+z+WdD/wQmBjWxo1yLgFKknXAc9vcuhM4APAsc2KNUmLYdJHU2Y0dT0zsfirZC1AT0/P5HlF3Mw6LiLo7xvbiL6xkDQX2BYR2yXtR/G45IHxOt+4BaqhorSkFwGLgMbV1ALgVklLKa5g9illXwD8ItMXNEmnVGarpGnA7sC2TF82qMwNwKMUt/Sm5VVVs7qancfMrBoiiO3jP+pP0muBTwFzgW9Iuj0ijgNeDpwtqQ/YDrwtIraNVzsm/BlVRNwVEc+LiIURsZAiOCyJiF8C1wIrcyTfIooofXNEPAz8TtIR+YzpzcA1WeW1QGNE3xuA7+VzrPXAsZL2zEEUxwLr89j1mZcsW67rzTn67wjgiTy3mVmlTNCov6sjYkFEzIiIvTJIERFfiYjFEfHiiFgSEV9rS6eGUKn3qCJis6QrgbuBPuD0iGh80qcBnwd2Bb6VG8DFwGWSeimupFZmXdskfQS4JfOdXYr47wPWSToHuC3rAPgmsIJiEMeTwCnj0U8zszHxFEoTK6+qyr+fC5zbJN9G4KAm6X8EThqi7kuAS5qkP0AxZH1wegCnt9h0M7MOcaAyM7MKK5ai93pUZmZWVR0e9TfRHKjMzOom2vceVR04UJmZ1UzAhAxPrwoHKjOzuvGoPzMzqzYHKjMzq7JJNphCxatDNhaSHgF+NoYq5lBM7TQZuK/daTL1FcbW3xdExNyxnFzSt7MNrXg0IpaPnK26HKgqQNLGiOgZOWf9ua/daTL1FSZffzutoyv8mpmZjcSByszMKs2BqhrWjpyla7iv3Wky9RUmX387ys+ozMys0nxFZWZmleZAZWZmleZANUaSzpAUkuaU0lZL6pW0RdJxpfTDJN2Vxz6ZqxWTKxp/KdNvkrSwVOZkSffldnIpfVHmvS/L7pLpyrp7Jd0paUmb+nm+pHuzzqsl7dHN/R0rScvz8+iV9P5Ot6dM0j6Srpd0j6TNkt6V6bMkbcjPeEOujN0o07HvuE19nirpNklf7/a+dqWI8DbKDdiHYsn7nwFzMu1A4A5gBrAIuB+YmsduBo4ERLFC8fGZ/nbgwtxfCXwp92cBD+TPPXN/zzx2JbAy9y8ETsv9FVm3gCOAm9rU12OBabm/BljTzf0d42c1NT+H/YBd8vM5sNPtKrVvHrAk958L/CS/x48C78/091flO25Tn98DXAF8PX/v2r5249bxBtR5A74MvBj4KTsC1WpgdSnP+vzHPQ+4t5T+d8BF5Ty5P43ijXeV8+SxizJNmacROI4E1pfzlMpsAea1ud+vBS6fLP0dxefz5/Y1+4yqtgHXAK8sf3b5/W2pwnfchv4tAL4LHMWOQNWVfe3Wzbf+RknSCcBDEXHHoEPzgZ+Xft+aafNzf3D6gDIR0Qc8Acwepq7ZwOOZd8i6mhxrl1Mp/qIc7nzd1N+dVcU2NZW3qQ4FbgL2ioiHAfLn8zJbp7/jsfo48F6gvCRut/a1K3lS2mFIug54fpNDZwIfoLgd9oxiTdJimPTRlBlNXSMarr8RcU3mORPoAy4f4XyV7+84qmKbnkHSc4CvAO+OiN/mI5emWZukTeR3PGqSXg38OiJ+LGlZK0WGaEfl+9rNHKiGERHHNEuX9CKK+9d35H/uBcCtkpZS/HW0Tyn7AuAXmb6gSTqlMlslTQN2B7Zl+rJBZW6guG2wh6Rp+VdZs7qanWdU/W3Ih8GvBo6OvGcxzPkq399xVMU2DSBpOkWQujwivprJv5I0LyIeljQP+HWmd/o7Hou/Bk6QtAKYCewm6Ytd2tfu1el7j92wMfAZ1WIGPox9gB0PY2+heODfeBi7ItNPZ+DD2CtzfxbwIMWD2D1zf1Yeu4qBD2PfnvuvYuDggpvb1MflwN3A3EHpXdnfMX5W0/JzWMSOwRSLO92uUvsEXAp8fFD6+QwcYPDRKnzHbez3MnY8o+rqvnbb1vEGdMNGKVDl72dSjBbaQo4MyvQeYFMe+zQ7ZgaZmf9weylGFu1XKnNqpvcCp5TS98u8vVl2RqYL+M88x11AT5v62Etxv/323C7s5v624fNaQTGa7n6KW6cdb1OpbS+luNV0Z+n7XEHx7OS7wH35c1YVvuM29nsZOwJVV/e12zZPoWRmZpXmUX9mZlZpDlRmZlZpDlRmZlZpDlRmZlZpDlRmZlZpDlRmZlZpDlRmTUhaKOkPkm5vU32H5OwIbZHLdPyfpJ521WlWVQ5UZkO7PyIOaVNdh1C8VNuynI6nqYj4G2DjWBtlVgcOVDbpSHpJLrI4U9Kzc/HAg0Yos1DFwpGfk7RJ0uWSjpH0w1z8bmnme7akSyTdkgv1nZiL4p0NvEnS7ZLe1Cxfln+LpKskfQ34jqR5kn6Q5TZJetm4f0BmFeNJaW3SiYhbJF0LnAPsCnwxIja1UHR/4CRgFcW8b39PMR3RCRSz6b+GYvqd70XEqSpWQb4ZuA74MMX0Tu8AkHTe4Hw5ez0U6xMdHBHbJP0zxVpF50qaCjyrHZ+BWZ04UNlkdTZFsPkj8M4WyzwYEXcBSNoMfDciQtJdwMLMcyzFbN1n5O8zgX2b1DVcvg0RsS33bwEuydnO/zsi2vLMzKxOfOvPJqtZwHMolmKf2WKZp0r7/aXf+9nxR5+A10fEIbntGxH3NKlruHy/b2SKiB8ALwceAi6T9OYW22rWNRyobLJaC3yIYgHINW2sdz3wT8qFyiQdmum/owiKI+UbQNILKBb++y/gYmBJG9tqVgsOVDbp5FVJX0RcAfw78BJJR7Wp+o8A04E7JW3K3wGuBw5sDKYYJt9gy4DbJd0GvB74RJvaaVYbXubDrAlJCynWLhp2NGAnSboBOCMiPEzdupqvqMya2w7s3q4XfttN0vUUi+893em2mI03X1GZmVml+YrKzMwqzYHKzMwqzYHKzMwqzYHKzMwq7f8BvE81qTHTdO4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds.sel(time='2018-11-05').u.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will generate some test points for interpolation.  These need to be in projected coordinates, so we use the `transform_coords` function to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0, 38758.62082587721, 57921.845123659696] [-555818.1242497492, -443013.06322267157, -328491.107208524]\n"
     ]
    }
   ],
   "source": [
    "lat = [85.,86.,87.]\n",
    "lon = [0.,5.,10.]\n",
    "x, y = transform_coord(4326, 3408, lon, lat)\n",
    "print (x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot as a sanity check and to make sure the points fall within the region with data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAGKCAYAAADQTkKUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOy9eZwlWVnn/X3OibhLZtbWKy3QdCsIiB9A7QEdxxFFsd1xm8EFedWRcXv11XldUEd4VV5Rx413RGUEaRBEBH1FlKVFEXdpFkE2G2iWll6pLbd7b0ScZ/4450TEzcyqyqrKyqzKer4fLvfeuBHnnojqvL94nvMsoqoYhmEYhrE3uL2egGEYhmFczpgQG4ZhGMYeYkJsGIZhGHuICbFhGIZh7CEmxIZhGIaxhxR7PQHDMAzj8uHBMtYJYUfGup/Z61X15h0ZbA8xITYMwzB2jQmBr+O6HRnrt/jIVTsy0B5jQmwYhmHsGgJ42aHB9kkZDFsjNgzDMIw9xCxiwzAMY9eIFvEOmcT7xCI2ITYMwzB2lR1zTe8TzDVtGIZhGHuIWcSGYRjGrrGjrul9ggmxYRiGsXuIuaY3Yq5pwzAMw9hDzCI2DMMwdg1zTW/GhNgwDMPYNXa0oMc+wVzThmEYhrGHmEVsGIZh7CJirukNmBAbhmEYu4ZgrtiN2PUwDMMwjD3ELGLDMAxjVzHX9DwmxIZhGMauIVbQYxPmmjYMwzCMPcQsYmNXEZHrgfcAh1S12ev5GIax+5hreh6ziI0Lioh8WES+KL9X1Y+q6tLFKsIiMhSRF4jIR0RkWUTeLiJfumGfJ4rI+0RkTUT+UkQe0vvsC9K2EyLy4S3Gf6yI/HX6/E4R+akzzOdnRORdIlKLyLM2fHadiLxaRD4uIioiN2zj/L4pnduqiPz/InJF77MXichMRFZ6D3+KcQYi8sr076si8oQNn/+wiPxLuoZ3iMgPn2FeN6Trtpau7Rdt+PyU895irKGIvFBETorI3SLyQxs+f6yIvDV911tF5LGnm5uxs+SCHjvx2C+YEBvGPAXwMeDzgUPAfwdekUVORK4C/jBtvwK4Dfj93vGrwAuBUwnPy4A3p2M/H/huEfmq08znA8CPAH+6xWcBeB3wdWc+LRCRRwG/BTwVuBZYA563YbdfSDdKS9u4Yfob4FuAu7f6OuBbgSPAzcD3ichTTjPW7wFvB64EfgJ4pYhcfRbz7vMs4GHAQ4AvAH5ERG5OYw2APwZ+N83tFuCP03bD2BNMiI0Lhoi8BLge+JNkXf1IsnxURIq0z5tE5GdF5O/SPn8iIleKyEuTRfOWvqUnIo8QkVtF5KiIvF9E/tNOzllVV1X1War6YVUNqvoa4A7gs9IuXwu8W1X/QFUnxB/9x4jII9Lx/6SqLwE+dIqvuAF4qao2qvpBopg96jTzuUVVXwssb/HZPar6POAt2zy9bwb+RFXfrKorxJuJrxWRA9s8vv/dM1X9VVX9G2CTWKvqL6jq21S1VtX3E8Xvc7caS0Q+FfhM4Jmquq6qrwLeRXeDcbbz/lbgZ1T1mKq+F/hfwP+RPnsC8WbrV1V1qqrPJd40fOHZXgPj3Mi1pnfisV8wITYuGKr6VOCjwFcm6+oXTrHrU4jWzgOBTwH+HvgdotX4XuCZACKyCNxKtCqvAb4ReF6ymDYhIs8TkeOneLxzO+cgItcCnwq8O216FPDPvXNcBT7IacR0A78KfKuIlCLycOBzgD/f5rHny8a5fxCYEc8v8z3pJuetIrItS/tMiIgAn0d3DRGR14jIj/Xm9SFV7d9s/DPdNT3tvEXkx0TkNen1EeCT+vtvMdY7VVV7n7+T7f/7GTvAbrmm0xLFvSLyL71tzxKRfxORd6THl13Ic90OJsTGxcDvqOoHVfUE8Frgg6r656paA38AfEba7yuAD6vq7yRL623Aq4Cv32pQVf0eVT18isejzzQpESmBlwK3qOr70uYl4MSGXU8A27UqX5Pmuw68D3iBqm7Xoj1fzjT35xJdutcQrc4XiciWVuxZ8izib83v5A2q+hWq+pxtzuu0n6vqc1T1K3r7smH/bY9l7DteRFwa2civqOpj0+PPdnlOmzAhNi4G7um9Xt/iff5xfQjw+L5lS3RbPmCnJyQiDngJ0fL6vt5HK8DBDbsfZAvX8RZjXkFc0/1pYAQ8GPgSEfme9Pm7e0FSn3ee8/+83ljZEj3t3JMr+RPpJufPiDchX3ue8/g+oqv4y1V1eordznRNz+aar/Q+P9+xjAtAzCPeHde0qr4ZOHrhz+r8MCE2LjR65l22zceAv9pg2S6p6ndvtbOI/OaGCOCVLcRpq+MEeAExMOjrVLXqffxu4DG9fReJ7vRTjtfjk4FGVV+cxO5O4OXAlwGo6qN6QVJ/vY3xTomq/nVvrOx23Tj3TwaGwL+eahjikt45ISLfDvwY8MR0rqfi3cAnb1jzfQzdNd32vFX1GHBXf/8txnp0+jfOPJrt/fsZO8RFEDX9fSLyzuS6PrJDp3XOmBAbF5p7iAK0E7wG+FQReWpaYy1F5N+JyCO32llVv2tDBPDSFuK0Fb8BPJK4tr2+4bM/Aj5dRL5OREbATxHXHN8H0ZJO28v4Vka9iNx/Tdu+Ke33AOA/M7+eOUc6xxHxb7VI4/ne5yOiKAEM0/tT8VLgK5O1vEi0zP8wr82KyNeLyFKa25OIEdGvPs3c+t83SHOT9Nk3A/8v8MWqeqrANQBU9V+BdwDPTGN8DVEcX7WdeW/Bi4GfFJEjKYjuO4kuSoA3EYPLvj/NP3s7/uJ0czQuWq4Skdt6j6dv45jfIN48P5Z40/ZLF3SG28CE2LjQ/BzxR/G4iPzf5zNQ+uF9EjG46+PEtJmfpxOi80ZiTvB/Jf6R3t2zoL85zeE+YjTvs4FjwOPTfDL/kehO/zNixPg68IZ07Emiq/cH07HvAP4ljXUq/lca4xuJaT3rxMC2zDqdO/Z96f2WqOq7ge8iCtu9xHXR7+nt8gPAvwHHgV8EvlNV33Saub0/fd8Dgden1zmn+meJqUhv6V3D38wHishrReTHe2M9BbiJeF2eA3x9utZnnLeI/LiIvLY31jOJAXQfAf4K+EVVfV0aawY8meguPw58O/DktN3YBXY4j/h+Vb2p93j+mb4/ZRs0qhqIf1+Pu7BnfGZkPnjQMAzDMC4cNxZjfeahG3ZkrG87+r63qupNp9tHYvrja1T109P761T1rvT6B4HHq+rpctwvOFbi0jAMw9iXiMjvEXPHrxKRO4nekidIrKamwIeJHrA9xYTYMAzD2FV2qzylqn7jFptfsDvfvn1MiA3DMIxdI6cvGR0mxIZhGMaukYO1jA4T4h3gqquu0htuuGGvp2EYhnFBeOtb33q/ql691/PYr5gQ7wA33HADt912215PwzAM44IgIh/ZyfHMNT2PCbFhGIaxa5hrejNW0MMwDMMw9hCziA3DMIxdxVzT85gQG4ZhGLuGCDgT4jnMNW0YhmEYe4hZxIZhGMYuIohFa81hQmwYhmHsHgLOhHgOc00bhmEYxh5iFrFhGIaxawgg3mzAPibEhmEYxu4h2BrxBuy2xDAMwzD2ELOIDcMwjN1DxIK1NmBCbBiGYewq4swZ22fPr4aIeBF5u4i8Jr2/QkRuFZHb0/OR3r7PEJEPiMj7ReRLets/S0TelT57rkgs2yIiQxH5/bT9H0Xkht4xT0vfcbuIPK23/ca07+3p2MFuXAfDMAzj8mTPhRj4AeC9vfc/BrxRVR8GvDG9R0Q+DXgK8CjgZuB5IuLTMb8BPB14WHrcnLZ/B3BMVR8K/Arw82msK4BnAo8HHgc8syf4Pw/8Svr+Y2kMwzAMYweQlEe8E4/9wp4KsYg8CPhy4Ld7m78auCW9vgV4cm/7y1V1qqp3AB8AHici1wEHVfXvVVWBF284Jo/1SuCJyVr+EuBWVT2qqseAW4Gb02dfmPbd+P2GYRjGDiBeduSxX9hri/hXgR8BQm/btap6F0B6viZtfyDwsd5+d6ZtD0yvN26fO0ZVa+AEcOVpxroSOJ723TjWHCLydBG5TURuu++++7Z7voZhGIYxx54JsYh8BXCvqr51u4dssU1Ps/1cjjndWPMbVZ+vqjep6k1XX331VrsYhmEYGxFBvNuRx35hL6OmPxf4KhH5MmAEHBSR3wXuEZHrVPWu5Ha+N+1/J/Dg3vEPAj6etj9oi+39Y+4UkQI4BBxN25+w4Zg3AfcDh0WkSFZxfyzDMAzjPBGs1vRG9uyWQlWfoaoPUtUbiEFYf6Gq3wK8GshRzE8D/ji9fjXwlBQJfSMxKOufkvt6WUQ+O63xfuuGY/JYX5++Q4HXA08SkSMpSOtJwOvTZ3+Z9t34/YZhGIax41yMecTPAV4hIt8BfBT4BgBVfbeIvAJ4D1AD36uqTTrmu4EXAWPgtekB8ALgJSLyAaIl/JQ01lER+RngLWm/n1bVo+n1jwIvF5GfBd6exjAMwzB2AgFxZhH3uSiEWFXfRHQNo6qfAJ54iv2eDTx7i+23AZ++xfYJSci3+OyFwAu32P4hYkqTYRiGseMIbh+t7+4EdjUMwzAMYw+5KCxiwzAM4zLBui9twoTYMAzD2DXEhHgT5po2DMMwjD3ELGLDMAxjV7FgrXlMiA3DMIzdQ/ZXneidwG5LDMMwDGMPMYvYMAzD2DUEcFbQYw4TYsMwDGP3EPZVw4adwK6GYRiGYewhZhEbhmEYu4p1X5rHhNgwDMPYPSxqehPmmjYMwzCMPcQsYsMwDGPXEAvW2oQJsWEYhrGr2BrxPHZbYhiGYRh7iFnEhmEYxu4hIFbQYw4TYsMwDGPXEMSaPmzAroZhGIZh7CFmERuGYRi7h2B5xBswi9gwDMPYPVL60k48zvhVIi8UkXtF5F96264QkVtF5Pb0fOSCnu82MCE2DMMw9isvAm7esO3HgDeq6sOAN6b3e4oJsWEYhrGLCOLcjjzOhKq+GTi6YfNXA7ek17cAT97Z8zt7bI3YMAzD2DVE2Mmo6atE5Lbe++er6vPPcMy1qnoXgKreJSLX7NRkzhUTYsMwDONS5X5VvWmvJ3G+mBAbhmEYu4jsda3pe0TkumQNXwfcu5eTAVsjNgzDMHaTXYyaPgWvBp6WXj8N+OMdOa/zwITYMAzD2JeIyO8Bfw88XETuFJHvAJ4DfLGI3A58cXq/p5hr2jAMw9hFZFsRzzuBqn7jKT564q5MYJuYEBuGYRi7h4B4v9ezuKgwITYMwzB2Ddn7YK2LDrsahmEYhrGHmEVsGIZh7B4CbpfWiC8VTIgNwzCMXcVc0/PY1TAMwzCMPcQsYsMwDGP3EAvW2ogJsWEYhrFrCOxaHvGlgl0NwzAMw9hDzCI2DMMwdg9zTW/ChNgwDMPYPcSipjdiV8MwDMMw9hCziA3DMIxdxZlFPIcJsWEYhrFriOxe96VLBbsahmEYhrGHmEVsGIZh7CoWrDWPCbFhGIaxe1j60ibsahiGYRjGHmIWsWEYhrGrWLDWPCbEhmEYxq4hIjjv93oaFxUmxIZxEXLfyTUASgfeCQAHFsZ7OSXDMC4QJsSGcRFy9cGFVowzR5fXuOLAwh7NyDB2DgvWmseE2DAuUq4+uMAnltcIjbbb7j6xygMOLe7hrAzjPLFa05uwq2EYFzFXHligDkodlFmIgnx0ee0MRxmGcSlhFrFhXORckyzg4ytrBD3DzoZx0WMlLjdiQmwYlwiHl7p14/tOrnH1QVsvNi49xFzTmzAhNoxLiH4Q113HVwFokpkcgOuvWNqrqRmGcY6YEBvGJcaJacPBgaMJ2qY2hfTZB+5bxgmtC/uhVx8A4I77lwG48aoDuz1dw5jHSlxuwoTYMC4xHnr1Ad5/70nGhYOgVEl16xDdfn0+lARYeu8FUKJYa2/N+WHXmEgbu4OtEc9jQmwYlyAPv+Ygd9y/TJPe9zKcaJJ53Dc6svC6DULd5333nEQEXJJtE2bD2B1MiA3jEmV5FlgoHU5g1lNiJ+BFNkVYq3aCLdJZwwELxTZ2ERHEWYnLPnvmHxCRB4vIX4rIe0Xk3SLyA2n7FSJyq4jcnp6P9I55hoh8QETeLyJf0tv+WSLyrvTZc0Wig05EhiLy+2n7P4rIDb1jnpa+43YReVpv+41p39vTsYPduB6GcbY8+pMOxfziJuYZN5qeA1Rpe6Pd51V6NBofgfQ6RCs6aGdNG8YFxfmdeewT9tJRXwP/TVUfCXw28L0i8mnAjwFvVNWHAW9M70mfPQV4FHAz8DwRyf8SvwE8HXhYetyctn8HcExVHwr8CvDzaawrgGcCjwceBzyzJ/g/D/xK+v5jaQzDuCh5xLUHecS1B3n0Jx3Ci1A1ShVCK8xNoBXhtjBI04lvE6BRTdZytIyzUBuGsTvsmRCr6l2q+rb0ehl4L/BA4KuBW9JutwBPTq+/Gni5qk5V9Q7gA8DjROQ64KCq/r2qKvDiDcfksV4JPDFZy18C3KqqR1X1GHArcHP67AvTvhu/3zAuCn7oj/9l07Zn/Ol7+KsPH+Xe1Sn3r1Xcvzbj6FrF8UnFWtWwPK05tl5xYlKzVjWsV4H1KrBWNazOGpZndSvGQbHCIcYFRMC5nXnsEy6KNeLkMv4M4B+Ba1X1LohiLSLXpN0eCPxD77A707Yqvd64PR/zsTRWLSIngCv72zcccyVwXFXrLcbaOOenE61wrr/++rM6X8M4V77keX/Lox54iC/7zb8D4M++69/z5N/+Bx553UGOjEuWZw3QMCocToTSCaXvIrSipesI2uBECKqU6QetajTue5qALsM4bwTE2iDOsedCLCJLwKuA/0tVT8rG/Iverlts09NsP5djTjfW/EbV5wPPB7jpppvMfjAuKJ/3i39JaJSlg0Pe+bHjDIoonq9+z90AfOi+Fa6/MlbaKp0wLTylz0IcRdkLSXwDLsSALoCqaSi9MPSeoIqqKbFh7CZ7KsQiUhJF+KWq+odp8z0icl2yhq8D7k3b7wQe3Dv8QcDH0/YHbbG9f8ydIlIAh4CjafsTNhzzJuB+4LCIFMkq7o9lGHvCv/t/3sBooQRg5eSUn/qaT6cKyuvecw8Ax5enfMvn3sA9y1O8EwaFY9AEfLKIF0pP4R2lk1gApIli7SS+dwJD6dx8IjH16QP3LbcFQQxj55B9FWi1E+xl1LQALwDeq6q/3Pvo1UCOYn4a8Me97U9JkdA3EoOy/im5sZdF5LPTmN+64Zg81tcDf5HWkV8PPElEjqQgrScBr0+f/WXad+P3G8aesXZyytrylJ988qM4Ma1ZmdX8h4deyYlJzVfd9CDuWZ5yYq3ixFrFfSen3H18wtHVGcfWK5ZnDSuzmpVZw8o0rhEvzxrWqqYrj5mCs7wTvOus5TuPrrQPw9gRBIua3sBeWsSfCzwVeJeIvCNt+3HgOcArROQ7gI8C3wCgqu8WkVcA7yFGXH+vquZ6Bt8NvAgYA69ND4hC/xIR+QDREn5KGuuoiPwM8Ja030+r6tH0+keBl4vIzwJvT2MYRst0dXnu/XDxwlqNayemTCcVv/4D/4FpE1ivGqom4JywVjXcd3LKtA7M6ph7NEs5SAPvWBj4JK7xUToXXdTpdekd3klMW/KdAOdc5D65tjXAdYetJ7Jh7BR7JsSq+jecOizkiac45tnAs7fYfhvw6Vtsn5CEfIvPXgi8cIvtHyKmNBmXMdM3vZTimgfC+CA6GKPlGC1HaDkG2V1H0rt/+Sv56Vvfn6KdY4QzwMpqzcqk5t7laXw/qQBSLrGyMPDMmpJB4RgPPIPCMXMhri83MPLKtHYEVZzAtIku6VJddJU5RURw0Na0NozzRfZRG8SUCnsmgqoeP90Oex6sZRh7zfKLn4Ur459CqKLIDa5/GGH1JC40uFCjoY6BTH4AoqDR6pR6wmwWLcXBkQdcsDk+56d+nef+2n/j5LTm6MoMgLVZw4n1ijuPrtEEZbIahTjfJ5woPeNRhXfCeOYZFq5bQy4c605YrRpGhWdUOMalp3TCqHBcd2DIwDsGPq4nkxpMZCP5zqMrPMg6PRnnQnZN7w8+nh6nu1P1wGlTa0yIDWMLwspxqGdoXeGaBrcYcCGg5RBcEdUu1HPH1B9/P8UnPfyCzek177yLz3zIEdZmDYUTTqxXrEwqZrOG9ZUZGrRd63UihEZxTlifxRWcWR2t4SZo+3pWB1aoGSSRPjgsODgqKdcqjoxKwKX61EIIihdBrdiHYWTeq6qfcbodROTtZxrEhNi4rDn+/B/HjwZMPnECDdHKLUZD3PIxdLqOqyuoKwgBt9jgyhFaDEEcKrGogFSTtnBzfdftFNc9bMfnufb3z+U5f3k7J9ai1bs2a5jVDXd/Yo2mVqppTdOrT+lTx4e6dCxPKpoQGBQ+lrysQ7tmXPRczoPCsT5rYhWuNFbQAlXPwEe3tcrmDk+GcXbsq6jpz9mJfUyIjcuKEy/4SarVddbvjUs2B66/llDVTI+vEKq6dVFrE3CDCcVsgtYV2jQQGqQcI80suqj9ABpaN3X2CVf3fZTy6utp7ngbAP7Gz9yRud/6rrt55AMPxfNYnzGrAyFEEa6rQFN3Qty4dGNQ+dYiboLSBD+33ptf597GAx/FeOGaJU6kNWeAgGOAgIueecM4H/bLGnGKQ0JEPgW4U1WnIvIE4NHAi1X1eN7ndJgQG/uW2d++or3zHnzO17H2qv+xaZ/Vuz5BqGrqSVx39WWBNoFiPMCPYr+PslhFByPCZBU/XUXCEFyFuEkss9c04D0qDh3GCOrpyaPtH9dkfZ3ReHxe5/JF//NvmKxVc9uOr1WERpkmwcwWsUsma9NEcc4W8FY0QZn2BHxh4FmvGsYDz3WHRjiR5JouUHWoh4E3k9gwNvAq4CYReSgx0+bVwMuAL9vOwSbExmXB5HXPB2gtXlcWzJbXmN17DOcdkh7aNDRVjRsUuCagTUDrWVwvrip0fRVC0wq8lMNoEddAMUL0BM3S1QCsX/tIiqRZH/nECg+5cok77o+pTzdetTnl6cFPfREPf/xDWRoVPObBh4EotnceWwPgpodfzcq0TpatMigcRenw3tEQWgEGECeICBoUDZ07ui+6TVBWJjVNHdp132XvOLBQcuXigPHAU7oYVT2tA4dGBQulR5Llf/u9y3gHn7zFuRjGKZF95ZrOhFRG+WuAX1XV/287a8MZE2Jj3yKDEYQopG7xIDqbRGFNluP0+DLV6ox6vaaZNYSgLF2zyPBIlyMr3iHOoyEgoSGsnURmE6QokeGIsL6KjBeRQXRZN+OrQRyS3NUnZ4H1OopcvyjGq99zNycmNR89tsZPPPFTWfwPP8gVNz6m/fxvb7+fpVHBQ66Mc7nu8Ij7liccWohWunfCJ1ZmNLUyWhzQ1AHva0KIAVoARekpBr7dH2gLeDRBWZ81rByf0DSBkAS6KD2zac2dA8/arGHlYM0VSwMODUsmdeDAsODapQEDHwt/QHeTYRjbY3eFWEQ+DCwTF5JqVb3pAnxNJSLfSCwC9ZVpW7ndg02IjX1J/c9viAIKaFVBPYPQUCyMaSYzQlUzOLCALwuaxZqQ1lH9qKRanRBmNcE7wqyGIv09OY8UA2QwQmcTKMoowkUJoUaLEYigroh5xwqNwrFJzcA7Dg4cdYBffvOH+NJPu7ad6+Azvp1y8RDX3PAA1ldmHLtnhYWDQ6ppzfosuolndeCKpeH8OdaBcuhRVZwXxEULGKJF7JxQlI4ipSv1aZKlrEEJaa0ZYoWt0Chrs4alUWB91rA+axh5xzBFXKvG2LSgKbUJS2cyLnq+QFXvv4DjfxvwXcCzVfWOVP3xd7d7sAmxsX9xPrqRM8UAAfxoHV9VDA8fIFQ12gT8aEAxHlCtTqhWJxTjAcPDB+I6cQgxuKSeAYtIWYJzuKXDaDFEU+BJGB0ijA+1EdRNUERiz9+PLK9z/1rF3csTvAi//qYPcuLoOsfuiVbydY/+PKbrNSePrrN4aMhsvY7W7Dhas0ujeDPgRTgwKrh3ecq1h0dM60DhhPtPTvFFaIUYohh7H4t5wGarGKLw9qOt+8fn9eNZHahSalTopS7lZee8ZHzX8VWruGWcmX3YfUlV3wN8f+/9HcQqkdvChNjYn9QVmkRYvEcZIIPopvajAUPvKEbDNiArrg+nlJ1U1MOPBvjFA2jTIAsHoZ4ho4Uo6EsHCGW0UGORD0ezeCU1jjpEU7gKUDedcHmB8cBzYFRweKFkbXXGYFzykH//lfjCsXz0JOMDUciKgWfp4JDxwHM4uaOBNg/YO2E8KBgPOqH+yP2rm4V4Q5pSvUWj4dzxzLloVfcp2tKYkgK3BJ8eks5J2oAuw9gOspO9hK8Skdt675+fOuP1UeANIqLAb23x+Z5jQmzsO8IH/mHeGs6RzcR1Y+c8Us0orhi326hnaAjodJ1i2EU4S1HiBiN0soqMFiGEdlxJ7mgdjGMqE9EVDTBplLpR1mvFu7i9SiK4NCo4tDDg/uGUumoYLQ5YX5kyWlzAe0c5LGIQlpO2VjTAuPQsDDzLk5px6Zk1oQ3cArjm4JDxoGB9VrM2axgWjuVJ3YopdEKcj/He4YbSWsXexwCwhYFnXPq2POZC6TmQnr0jPaQVYZE97CBjXM7cv401389V1Y+n3va3isj7VPXNuzG57WJCbOwrmjvehkB0H0OMeE7IYBTXdUODhoAbL3Z35km4tZpPEZKyRIoBAXDjRdziAXR8EPWD+BiM49qw8xybgZeAiFA3SqMQUJanUbirnnA2QRkOPKOFkmraMByVMdI5uZNH45LDCyXjQfwTXRoWbVOHHDE9dnHtmOTly+I8KDq333WHx6zPugpghZM517RPa8dSS/veecegiKJ/aKFksYyiPCwcC6VnmPoct9covbzm0Knd0tYwwmjZ5RKXqvrx9HyviPwRsZfAjgixiDwDeJ2qbjtCeitMiI39hyugGEDTRPGFnoUcBVqIwiyDEVJ2rl+tZrGSVv6hcC4K8KErEediEFYxigFZw0XUD2iKEVVQZk2TXLXR+lXi+vC0CW23pPVZw8okub6dMByXbZQzQDksGC2WXLk0YDwoouCWnhPrFePSc3ih7Los1YEmCWuumBXLVzat+xrgioVVrTMAACAASURBVMUhR1djY4i6dyMQglKULkZapzU754TxqODQuORQ+q7/8riHzF3ejx5dacc+lSV894nVvFS+qRKXrSVf7giyS0IsIouAU9Xl9PpJwE/v4FfcAfyAiDwG+Gdi5783qOqxsxnEhNjYN1R3fxBXlGhdwXARGcYfe6dhvi60OHQwJpQLsQhHz40toW4rZWl2N0MXDT1YAF8SfBkDmRqlmjQ0qlRBWW1iVLF3Maq4Djq3fjooHEujonUR94ttDIoYmbw0KnnQkXErmNM68JArF7hiaUDpHMfWZm05yvVZw6wOrM2azoU9KJjVDd7FaOkDowLvJPYink2Z1TGoS1VjilOKqvZOOLww4IrFAZ923QGWhlv/PHiRs3ZFZzHO4pxTuSzS2rjAXAv8UYqDKICXqerrdmpwVX058HIAEfkM4GbgD0XEA39OtJb/6UzjmBAb+4qYQpQkQlwUXQ1dGcq0Xf0AHR+K9aIh5v1qQGFTyUp1RbtPjJL2TOvQrgNXQcl1MuqgFE5aEfYiVE2MNi6dYzG5mpu037iMlaxmddO6lK9cHLA0KtqmDN4Jh8clC6VnrWo4NC6ZNYGlUdGmNg0mNcuTmlkTYjBVKleZBTZ3XfJzLmVhmNojLgw840G0hK9cGnBwVHIgRVu/8+MnePQnHWqPc0mE5TRFp7PgeiG69YnPEG9QMpaDfJmySyUuU1vbx5xxx535rrcTe9j/nIgcBL4Y+C+ACbFxGeHLKKTFMFq4Ipt6B2sxaKOcQzmiCRo7FSXl2FgKMubLRjdzKRAQZj0RngUlZ/9MmhD31WgxRhd1aMcepvXYoMrgoGNlUjNOlmx2Ky8l6/VgSle68cgCC6VneVrjU+TyyIe4/qxK1QSqMh43WIn1p2d14EC6ecjnU7gY6Xx4YRCLeaTzO7xQMig84xScdeVStIgPDDyHRgVL6cbh9nuXedg1sYJW361874lu7TeT14P7Oh3oxDlX8Wo2B3AblwOye67pvUJVTxLLXr5qO/ubEBv7Bi1HnQu6GHQWcXpGJFrCItRBWZ82rRgEjZZeTseBaNGKRFcsQPCC6rwl3ARoVAka93NOaFQpUorPJKl06R1V1VB64cg4urVL57hqsRdpXUdRXihjacmlgWcplZkcLgwIqgx9YNrMl6nMOb5eJFnXoQ3oyq/XZl0+9eGFAQuDJhXtKNto7BjNXbI0KFJQlqdwXWDWHfcvt6+zS/lUAVpZhPsCnK9z3m5cruzLEpfnhQmxsW/QYghka9iB8+0acBOUOiihjoFNqspKFYgpv4qjCz7yQmsleyc44nNryYUYjJWt4X5Qkgh4pF0DrkIshjGpA6UTylQXcugd4zIKqBOhbkJMDZIYMHXN4oBRkff17VilFxaCpwqdlAWNEdkAC6VnUgeqEC3jpVHByqRmIbmwx6VnZVozKHzrvo5u6fgYecdCGdeqy1TG0rvTdz0/5b/HBos3aGcNZ08DwIfuX7Z61cZljQmxsX8QB00VXdSuIBDFdFaFnmh2a7rrdaBqUp5viILjUiCST/mxPm2XkEV5PiK6SZZoo4onWo9VUFBYqxomySIFGBY+rRXH/VI8N42CH3QWQqOd5Tj0ntLHOVG7dDyMUqhUbvQwTXnSa66hdMKkEUbeMWnC3BoxwBIF61WTzllScZCcouRbixyiYDYBCrfZbX86VDsX9r8dW20FOFvIWYTNO30ZIuzaGvFuIiKPBm6gp6uq+ofbOdaE2Ng3DBcPsLIWVz+ndbQ262SNNgqzRpk1oRXItSoVsUj6UnrXRjgfGBSxilRPgPoC0sz5VpVB0VnBqrQBWhvpW8Xd+1x6kraRghPhwKCINwECDmHk4+tZv1qXg4ETwHNgWDCpAsuzmmkdWJ41DOsYvLU0KrhiadCmNcVqWY57TsZWqblox7Bw7XlDvOkIAvHXU9HTBGhl8k3ERz6xgsj8dYtXK46WT+OhV5s1fDkhyL4rcSkiLyT2IH433cqLAibExuXH0sKYE6vreInrp5M6ME3rubNaWasaqhCY1qF154bW0outBEsvVCG3FYwlJYNqWgvWTYJZOscsdVjyIgSilV2qgwKq0rcuaO+k/d40QvudC8M46KQOjArXuroHPUvUpzVoTWvS+UbBi7BItuSFFalboatCzBXGp7Xq9P3DwvOoBxygCsrQOw4MC9aqaFk3AXDJZS+xMEmj8bvPxMZ94rWLrzcengPADOMS57NV9dPO9WATYmPfEd3PsF5HF3QW4WkTRXitipHMVdCYXtNo6+ItHVQNjNJfRo56rhptBSYoKbk4BniNinj8qHA0mgK8UoQzNZQ+UDrfWoVOXGsxO4miXjrXrksvlD7NRSgcc1WsIAWPSfzuLNKxLaFj4DTlMPt2Tvl7c4clJ3EtuvCOqxYGrFVNG9F91UJMkxoVjtVZwPlcFlNQ0dbCPR1zkejxf9EqJt5AZKP64dccPJd/XuNSZ5cra+0Sfy8in5aaP5w1JsTGvuPw0gIQo3zzj38V4nowpNzW1MKvFeO0xht3iZZ0ti47aziOny3k7GKe1DHQKluRXqS1JJ0Ih4YxFSmKf8AFwUscrPRdlHRAKcXNpf3ENWuhdJ0rN9cBc2RrON5AeJU2wtuNYsT2pAptZHXVxMCxscs3DTG4K7vGS+daT0BQGJcu1bOmtegz//CRowB89kOu2HT9FeZuWrRXZ7vdwbiM2ZdR07cQxfhuYEpay1HVR2/nYBNiY99SOqEOUXh9cv+6IFBAmYKp/Ibay5mQBHkrV6x3sUpFE5Sqb8HmoCrSuq5EMYPOIhxqtIazpT30vmeNS7senN3bhctNFYQijV04mUutGvj4aBs5iMM7ZRyE9d7NQo60rhql9NJa4QeGRWvlDnv1suoU4LbSq1Wdz7MvyhtpAnOWc3alp3IpbYS6YewjXgg8FXgX55CdZ0Js7Fv6fw1RPDzBKaVKL+hKmdJPBUpR0EFpYG5tF7oevC5HeCVRad3OSCvCsV1gN4d+NLR3jkkqxyUSXcs5dQpAkdQ1qRsnC7IXmDW07Ql9qGBW4Xy0vAvvWCgclcKgClTBMypdGzSVl6hLL63Qh55oOqRdD88u/GxNAyz1Irxf+757Whc4wOd/ylVx/kl4c3pX/4amMZP4skf2X9T0R1X11ed6sAmxsW+5/oolbr93OQmdQ13fZRrXfaOwurmAoiYV99gY9RzXg7sfkOyyzRZiFvEswiLRys3jiAqaxq5CHKt0joHv1oJFYr6yJEvXSa6K1aVWiSqFV9z0ZCzHmWpiU8fGDgIQaobiKIcHYl/kkNO2euvcdOvQCkzqOP+Y/hRFdFKHFGHeUDdd96d8E7CRt37s+JwXoB/glsmu8Ld+7Dif9eDD5/Ava1zSyL50Tb9PRF4G/AnRNQ1Y+pJhADEq9477l3tbUioSUZzrEAV52jStsHbWW2/NNT0vlH6uJGZfYAbeRQu1Z8VKsqbzkLmoRc5FHjhJ1nDcP9aqdm3UtGuqmBvdP4NmBhqQKv29z1LBSl9CaNrP0YBbO8bVg0V0fIjgy/Z8m41WagAvcV5OQJM32rtYLGSaXNshaBtZDTnlSxgWjqsWyjYgLudWdzc80o4X1J3WtW0YlyBjogA/qbfN0pcMI3PjVQfaohK5AQFA8NndC+DBz9eGzvSF1/WEOVbeivtkQU4BzXFdOlmyg2QBVq1oC0VaNx14iUGkyfJdKCR2gKpifm/bDapJyugc0tSxbGfqGiUaYmOK0MR9Q43k/esKpwFtKsSXeF8wBMLoECCs1VF4V1PQmvjNa+Ybo66bFPw2q0NsJCFCFeJNyLRwc/nY0K1JexFiUnJoC4YYlymyv/79VfXbzud4E2LjsuCBR7qayHenRgVliuiV5CWLetM1Zuhbbfl9v52fS52VvMtr0KQI5mgF5/XXQVLrgYdZrj1NKqohSWxDE1OinIemmhPg+Dr7zeOT5BcaoihrSOMkEQ5NfC8u3nFUa0gYoHWsu+3rGTpYYDxcYtYooyKuO4c24EsITnEh1sZeq6ZtDvLqrG5bNK6nphVZkNer2FGq9G5ufbxqIDhwonjntpWPbOxXZN8JsYjcAvyAqh5P748Av6Sq376d402IjcuOMrUp9KIUDuoQLdcGbdcvs1HYWr7EdeAc8dugrUXsUoOHnO4E3bpu6VMwVSq9OSxKJK8ZV+tIqJHZWhRaDV1zimIAISBN1bmaoW1kMfc+vZZmBnUVWzqGGq0rpCiRWQBXxM+8j32V/QCZruBcQVmMUO0s4qlC6AWwrVWxxnX2SDdB22YSsxRwlstnHhwWNArDQts15LieHPO1S+/aNC+ILRb7JUUBHnGt5RcblxyPziIMoKrHUn/ibWFCbFx2+Ow+dpLkJhVdDPFVQLu14V4zCIekNU4oJOYcq3YpRuMiWrnZah74lHLUpHXbZoYrRnGwEFVNpqtIM0OaWXQvi0NdEa3avpXb76ec0pCiWEcXtOQ14yTAefxohkbRj67sAlzTBnlJPcH5Eu8chcY16oGPgVizWtvCI/m5UVpreFaHua5Oszqw6rr3Iecnq5AriNGENkDOn8Io6rdcNPYnus8sYsCJyBFVPQYgIldwFvpqQmxcdhxeWuDo8hpFbyGzSVWuorEmbYWsgM6JcF7TjSIsc7mxOZ+3cNl61mjx1nG9F1cg1VqMck7BVBLqJNJNrADgCiTUrbsZDVGEm9hfWeppK9hRpEMnuvUUrSu0mgEg5SCKssv7VOA8Ug6BVFejiSLvXVwjB6jqWLwkaAzMOjGpCapzXZ361nDRayoRS3DCtG6AWMErhpqF5K5O7nyJXoYcaQ5dkRVjnyPsO9c08EvA34nIK4l/Wv8JePZ2DzYhNi5LXC+YCvpN66Mo51SjbN2WKbp5aeBS+4P5tn7eCYtFFEqZTpA6imErktC2ZZR6BTSKtFs/0Vq1qgPEha4zTbaI6yhleUzRJP91+ryp0LoiTCdoPWubritEV3Q1gxDQdLwbLyKjuGYex5/iC8CVVMTuTk7idch1t3MRktznuB/QlUU4rxPnazqtmxi4ldbTXdC2yUUuKJLrZJc9T4I1gTAuNVT1xSJyG/CFxFuNrz2bcpcmxMZlSS6DeXR5LVXeAlJDBcJ8RHS/elXh5ot05HQoJ4LUk/SYIfV0zuqNUc4FIi6J8QQ3WZ6bUw7Q0pAc4iFWw5Bmtmn+ksQ77h8t3iy0mt3eNdFCDk3rzgbQponi7gcpSGwW5+bLJMDZOnUEbSido3SwVgEExgPfCnImFzyJ1y/euLi8Rixd16ncJCNX9eoaV5xbz2PjUkS6O9x9RBJeqzVtGGfLFQcW+MTyGj6buEgKJOrWhssUAV3knN8sdE3VFSZQoghXkyjCoUbqqnM/JwFGXGwBV0crFvpVhproOhbX5vHOiXBo0jpv3b7XEKCeRXFt9wtRfIve+7obR8IoPSe3ubi4Xpyt+2zRNg1N0LYhRFDFN9A42nKZG+kHrOXXhY950aV3DNPrPGZb/ISz63dsXOLsk/Q1EXmbqn7m+e5jQmxc9lx5YIG7T6ymGtJxfTQbe06iAJcuBXmFJgpwG0jl2sATqWdRLEONVNMuBQni/vkLW2O2QaeT6EIejtIXBsQli7gXoJVFu10TTjcD2jSd0OYgrvxZXSHOd++rCinL6L4Oo64oiCvQZoaEIYUrqEJcC18ofQrS8m15y1h1K+AlCmujSkhVu/L1itHiru1tnEV4oYyWcK7J7aXzMHgnbX73vx1bnUs3M4yLmEeKyDtP87kAh840iAmxYWzA9dY5pd/z1wmizVw0MyEgvojb6mmyijdEL4cGigFkq7QYdEI6S4FcoYGiTOlFPq5DJ2t3bv9q1h0Dm8S3HSsR15KT5Z2asUsqBKLZPa0h3kT4KcWgIFWoZOAdh0YFpRcmdUiBVJ7lFCnthVSHWik11qPuV9qKJTxlkwjndeGcd90XYdmHLktjHmVfRU0/Yhv7NGfawYTYMIAHHFrkruOx0EfuJxyjpSO56QIhB1DNupSiHERVT5DZOlpNO3dxslalqDp3XLWa3MpxXVermN8ro4VoHQ9G6HQSrdfpJEZMp3G0ruI2aF3V4j3aNEhZbj6x0LRuaylLpBhE0W4ahBlSlFD7aBWHWDykcAUjn9srxhaNQx+YFk3bqAJya8Mo0E0KyBpuIb6jZBnnHO2cEuZT+U9HJ8DZO33X8VWuO2xW8b5E9k9BD1X9yE6MY0JsGIkcsOXb9/PiEAOokus5V7vK/f40RMuzqdBpqv3cC5TS2SSWp0yuYp1NugCqJMpZqLWucKPFTnDTODqboFUVj83WtUsWtHOn7GmUhTqPJ4NRtI6L1CxCNVr3Lv4ceCeMS4dvostZBERi3m8U1pDqVUe3de533LeASydtr2PvOuHN5GYT+ec4N7oI2l1vw7hcOKUQi8hpF5cTlaq+awfnYxh7xjWHFrnv5BrQWcDZjZq7Hkk9ScU2UjRzv7JVimDOLuAczdwPTAmTVbRKruJs8SbXdbed1trNa7/ZctbZBJ2stpaueN+JcZHWkXudbcQ5NMR9ZDiKz8Wg2yeXx5yt40/ehfoBrhyztHCEpiiZNQERYdZoavygrNdR8vsNM4betxZvXv/NYgvRHZnrcQ+LTmkDsUtW9kaYFXyZsE8s4p3idBbxXwFv4fRZBTcCN+zkhAxjL8nro7kbUtt6MAVp5TrOQFtYA+hSlKATOddEN3OyfrOoivfodD1Zyxvcy77XHi5b0Pl9aKIlnKO2y8H894UNS1H97c6nYC3iPJyLBT+qKVJky7yGQSoqUk8pwiq+HKNukCKhHT7EutRBYSXVvCx90fZObt3OGzpQ1UER19Wxhvn63ybAlxP7xzWdEZHvA16aK2udLacT4reo6hee4cv/4ly+1DAuRo6vJGvY5Y5IqTpWv4lCv9SktKodc4PFIUXZpR4lcd2SFLylIUSx7e/nfE/Mfa9cZQ76mrd8o9WbrHHnNgVrSVG285nrAxuaVOwjCftwAZoGF1YI4lBfIi5+V+kHSCE0GvslNwreeeoiWsg5zcvn6mR06+xtlynp2kGa8Br7jAcAbxGRtwEvBF6vudrPNjilEJ9JhLe7j2FcSvRFWFJwFhqiG7rfBSkFnKhPdaEdaFF2aUqhQYsBwixGiWZxdK6LdE6CqjXIwCehdO1+4j2UMWJaszinRg6t6Obtad04NpDKIh6QAnCj7gRDA5TRQieeX1sSc5BSmkLf1V6nVCqldJ5CHOrjurAXmAGFl9jPuGfk9AOwvItFQgD0lCvZxuXEPoqaBkBVf1JE/juxH/G3Af9TRF4BvEBVP3im488YrCUinwu8Q1VXReRbgM8Efm2nosUM42Lh8NICy2vrrVRIv+4ltJ2RyMU5fPzz0SKlAKX8YgVEasTlStSr81/kHFS9Epi5JGXT9Ip7pK+uZj0R910P4izC2b3di46mqlohV1J+cw7MgrlcZA1NTGdKVjPTVRiMYxMKBkg9Q4uukpekphTOeQrfG9OTLN5OaLMQa67jHS8MVxxY2Na/h7GP2WdCDKCqKiJ3A3cTqwUcAV4pIreq6o+c7tjtRE3/BvAYEXkM8CPAC4AXA59/ftM2jIuPHLUbVBEnSBKPNigrWcFtxHOKNEYDWgyji9rN2rXcWN/ep7+0KHaxEYNPlm2KeJ5NorimFKQYZe278pR9d3PfBe18O16Org5VjXiHK8q0YjsvfBoaJM0l16UW77uiIb0iJASQatKdv2p0pfsCioax71Km1BXkkBJNbuigiorQBLVoaGPfIiLfDzwNuB/4beCHVbUSEQfcTtTOU7IdIa6T0n810RJ+gYg87XwnbhgXI4cWxwCsrK1H96vz0UUrDnEFmv9iykHnXpNknYqL7uu6RItR6qpUIotVtEpdbMAggxEMRp3IOodO1rr0pekkimpR9gp5VFHce6UqNfSqdYWGZjqlmcyoJ1O0CbiyYHBwkcHS4biPcyli2nX1qAfJbZ2LhjjfuuKlmcXz1tDWvUZDPA9foOUC6st23uIHKYBNWsuZ1M+5cGJOaSMi+7LW9FXERg9znmJVDSLyFWc6eDtCvCwizwC+BfiPIuLJt/aGsU/JoqGaoogBbVIlqj7Ot67rvlhLHV27ADJcRMphjFAejqNLGKJVmq3QwaitBy1lGYW3H0Wd60o73+Ul9+eZ3NL1JIpxtTrBjwb4QS9QawOSxmqFdDDqqoHlalw5vamZpZKZDophOvd0U9CQxHnQue5dgYrQV99999NrnDv7zzV940YRFpGXqOpTVfW9Zzp4O0L8n4FvAr5DVe8WkeuBXzy3uRrGpcEmFzXMie5c9HQvJzd6slOHJUdcZy1maJPaJ1bT6LIuB9E6Ho7S2rBvA7k0W77ZNd1GUScXdJh01bJyZHavCUO9HoV4enwZXxaMIY49l1+cA7pcl/ecbwqgy5FO5xqLicxSXewAQ9c2pFBXgErb1EJdQUAIwWxg47LhUf03yWD9rO0efFohToP9rqp+Ud6mqh8lrhEbxr5l43qmSsqQbbf7bs04u2+zOOdArnoWU5t8ip7WAilCl1KUKmIJed04Cd/MtcU9+uQ0pU35wlugTcP02CoHb3BtoFYWe3F+075S9LpApXXnnPakdRUjq7NQl4Mowv2mFMk1r+KiCG+RueH2nzvSOEf2S9R08hb/ODAWkZN5MzGh4PnbHee0QqyqjYisicghVT1xzrM1jEuA9clkruG90POsZms3iU/uL5xRumYK7baU2qRpP8nH+UEskzkYd9/VpifFCGcZjlJkc0/schSz8/Ni2nTVu4rRkLAYI5yHhw8wPHwg7tMX4TzPHKg1HMfAMT+I51dPU1OIDRcoR3f3052aAlxv/RhSje5OdBfGI9bWJ+1r43JH9k0bRFX9OeDnROTnVPUZ5zrOdlzTE+BdInIrvTwMVf3+c/1Sw7jYkd6zpsjpWN6D1mWrvQIV4nyb4kN/n2LQCbA4cLFOdSw72fsxCiXiim69ti7J+b6QqmE1TRttDXRCmvKKfV3hygLxjnJxjHiHH22xpu1c9z3DcbS0/SC5wrPbumeR91zphFhTO1v1sQtVV/1LQt2lefUwATb2IyLyCFV9H/AHW5WFVtW3bWec7Qjxn6aHYexrxqNOLFbW1ruiHglNYpxfb/S+tkFdga70pQbCMLmGc2UuDbEISHtgaGtXS1UixRCX835TFyepS+RQXE8OK8ejaKYOTa1YLx2G0ODrai66WnwUbxmMWuHNUd7qCtT5mAudo74BUkUtSaKtIcSSmqFBc5ZTKt0pGmJVsXIIqTzm4PA15/EvYexrYk7fXs9ip/gh4OnAL23xmQLbKnp1RiFW1VtEZAxcr6rvP6spGsYlzFZLmtoWqdhiv346U2/9OO+gveAnXK9cZhLmHGWt5SgKW5OCowogdUuSAmQ2iWu4w1EX5AW43NmprqLV2uRa1ummoByktemiC6zK71Pd7C5Kej4xQpzvuao3rF+n+UtTQ1NTXn39WVxl4/Jj/9SaVtWnp+cvOJ9xzng1ROQrgXcAr0vvHysirz6fLzWMix3XE1xVaIK2r3OhioxqrMI155qGueAtdUUX/ZxSffJjrmZ1OY77FjEXmdEBGC4io8Vo1aZcYClKpChxiweQwQh34AiyeBB/6Er8kWtwS4dx40Xc4sFoCQ9GnRUsLo1fQp6DL1MXp6Kblyu6B3SWdC+oa+4i5Cpfp2B9MmF9Mjm/fxjDuAgRke8VkcO990dE5Hu2e/x2XNPPAh4HvAlAVd8hIjee5TwN45Ija20W3Y2RwEE11qbuBSr1g5YAckOIOcTRtkbQEC3SlHucP48CXadiGoK6cRTAukKGVcw1bmJ3pzaga7jYBpSJc2jORw5N657OxUdyAFl7gwBdypUW7bpzPh9pxXi+xWIU5uKsCjT0xbi/HGBcRuwTi7jHd6rqr+c3qnpMRL4TeN52Dt5uZa0TMv+HZgmCxr4mqOJkcxpOLMQR6YuwaICm6taBs7sX2spcbZoTdGKsbKraJXVKdRKHNrPY+CG5ugWQ0VKq4JXWcVPEs5ap0Ebe13lwk66aV7aGnZuzwqFzm2s6bu5GgJSrLA4JdVwLTufQBaGd+YdVUrnLPitr6ywtjE9xhLFf2S/pSz2ciEjuuJRSfwdnOKZlO0L8LyLyTYAXkYcB3w/83TlN1TAuEbI4rKytb/m5E2mbQmRxbUU4B2LFih5dW6INKVAxujhV5QJI6UNapGjtVKu6FfZkKYs4pJrGfXKlL+fbql8aaqRIY7oqBlnlEpw9azjurL07iyIFjrku8CyfW9wBbXo/oDKfhlJee2pH2fpk0om+bvYuGMYlzuuBV4jIbxL/9L6LtJy7HbYjxP8n8BPAFHhZ+sKfOft5Gsb+QVKScbsu3CTrNLmXs+BqDn6iV46yF1Gdazvjyq70Zb8yV27om5+z4GpAi7KzSPM6NLTvBRAdR+s5B2i5aNXqxlKd7Yml1o79uWtI1bNCbK2YzrnddhbWzUYRtiIflyGyf4K1evwo8F+B7ybet76B2PxhW2xHiL9cVX+CKMYAiMg3AH9wdvM0jEuPpYUxK6k1oqTHluKx1Q9LFqlsdYqb60qEOHyvhJe0UdUpkhm6UpP9wK/BQhTI2VonnE0dA6VEUD+EJopxKBfmq36das6us9oRh5bpfROtfPVFO4+cerXd6OjJ+nqcV+7GlL8es4wvW/bZDVhq7vAC4G+I99zvV9Uzl8BLbEeIn8Fm0d1qm2HsSzaKMdD1Ks6IpGjnnqWZhS5HI6d2gNBZx3XQTtxzIJSmtWRxbenIvuUZfBk7Q40O4uppPCQ0SD3pvjPUaDGMIiplm6u8idzSMFUKm7PWQ51yiIkBYT6nPA0YHHnAWV3Dvgjn4mXWFtHYL4jIE4BbgA8TfyYeLCJPU9U3b+f4UwqxiHwp8GXAA0Xkub2PDhKbHu9bRORmQ4VKzAAAIABJREFU4NcAD/y2qj5nj6dk7DH9NWPp+Zk3BmD1Xb79cphZhBsFzWvLIq0YBdWuLKRIt96bLId+6c2mUVQ1irgMUFWKYkDpSyQJcxgf7kpOZlFNlmzfQlZfzAWW5apYmm4GpJ6lspxFay1vaV2fgslaKsYnPp1nd/6BeP4WsHW5sbuu6V36Pf8l4Em51oaIfCrwe2yz8cPpLOKPA7cBXwW8tbd9GfjBc5rqJUCKdvt14IuBO4G3iMirVfU9ezsz42Jgk1s6B2AR/08ltR3U0AY6IS5mAdGJEOl1o50gZzetF+ms5yS4GwlpvNzhqAqKFgVSFngnuKYi/u7EgCvREMU0vW/7IBddfjH01q/zOnC60dDWzX4KN/cGWgFO538q8tXMQXEmyJcHuxU1vYu/52W/4JWq/quIbLtd8CmFWFX/GfhnEXlZ2u9yqaz1OOADqvohABF5OfDVgAmxMVczebKeIqpzGg/RyosWc/dDUytkE1pEUNVWSL1EQQ4IPqlSEzrxbbTb1k8h7FvV8XiY1NHdOwAKXyKpuEacTpdKhUY3c7R0y00/itGCbuYFNwSQ0FnFidmxuwFO6apug74STuI1yq8N4wKzW7/nt6U14pek99/MvAF7WrazRnwz8D+If983ishjgZ9W1a8625leIjwQ+Fjv/Z3A4/doLsYlgM4FImnOSJrfh25ttE/ue6y9HNu836zJLtw0xhYDyIaU/qDxuKDgxVFIiuD2cZ1YNVm1vmz7BgOxnUW2nDedXy9POH0+XV1O30/7fktBz69TR6aYn83csWCW8GXH7rmmd+v3/LuB7yWm9wrwZrZZzAPOvbLWDWc3x0uKre7TN/0CisjTicW+uf56q617OTIad+KRg7lOh24Iltr4H1peBw4K06Z7nY9rtCs0Alnc0lgi4BSfPquDojEKjKJnsWdyz+BZE3ASj/fiKFJ6U86JbsW0CbExhCtQ33ncdLDYVemSrsBJkNiqsR/U1qZ89a9Jej6xuj53jQ4vLZzhahqXKpqWXnaIq0Tktt7756tqvw/wtn7PzxdVnQK/nB5nzblW1trP3Ak8uPf+QcT18jnSP/bzAW666SbLwTA2/cVvVZkLeuu7SuuO7m8XOqs4u53no47zh10wl6A4FVIm8rzbtxdxndeelWg1N+lROt2UUiJN7halXb7zXDOLTtwDQgjxJkGc3/RT1wTdMu2rCTqXwmUYZ8n9qnrTaT7f1u/5uSIi7+I0wq6qj97OOFZZazNvAR6W6mn/G/AU4Jv2dkrGxc7Swpi19VhDeWNtamXe2u0/Z0Rit+OgMAvaiyxmLuWpCV1tDxGiK5zs9o0WstKtOTcKhevWr9FuPtnazjfZmvbx0K0Pt6IfunXiXoESbVOu8nlp+9ohqRdzJ/6ZvqVfJevf9PgyQbfOpLtAXOjf86/YiUHOtrLW77HPK2upai3/u717j5errO89/vmtNbN3buRGwAAJBCG1onjBFOmptlQQonJAW9HYKlQ4RS3Wnnqpokc59dIWtdJjW60oWEEUqa2VejkI3tpTRbmoICIlghckSGJCSNjJ3jNr/c4fz7Nm1uw92Zkks2dmz/6+X695ZWbNWmueWdl7/9bze25mryZ8zxS4wt3v7HOxZBaZXBPOcqeWT96n9XU52LpPDcB56Y9XVq4Re7FmclgdCrMQfA1qeajp5iQknoGl5O6N8hS17cSd3EMATopgmTXXTo4lDh2voDlBSVTcFJTrBVnsoBZq897oIV4oauPl61G+OSlS1UsWqu14+HjPJnKZ6b/n7v6T4rmZHQWsdfcb49LBncRX6GRHdx8jBOK37G3fYeHuXwC+0O9yyOyyYP68lrmpizHDE5kXay2V/i2PIW4G26zoUe3lc0AMgTFF3fzMDKeaGG6xU1dM9ZZrugB1Eio4tUY62pu1ksmdp4q2YfdGL+vG+OgiOLthnjNvwUKguaJSEVzdPfYWb36/4jmENuzyxB6VeDFSm3qToqFNw6eXbXm9+HseV1q6AFgOHENIgf8DcEonx+81EJvZOuDNwJry/p3mvkXmkvKY3yIANYYg0doM264XtcfgXbQHF0G4CMAeg2hxngQLNeTcQjtxYo3zpjiZG7vreRhbbEZqTp2Q4nZvts860/xxLHf0aszclTG6eHljc1F7Tw121Zup9WKMdJY3b0YmMme0kjBezymWlKjnHlLpsR94Gkuzp0U3RAbMhYROzd8CcPd7zOzQTg/upOp8NfAG4A6aIylEpI3dmTd6Mhe/LJOndJzc8bFINbs3A1bRHlzUkBvvu1PLnDQJKfCMEIRzC0E18yIAhw/M89CTupiNK/Nm+2xRgy7ywpmBJUa1vPBDMeNWMRVmEmbYGj1oact3KNdWF5W2b35kjDym2j22f2cOOyfCGOei9l4Mo0ppdhqbQx1E55Qi+zNkxt19wpo/uxX2oeLfSSDe7O7X7WfhROaU7eNZc6pKoJKEYJMQglw6KbgUKdoiMBY1y1o+tRacxyBcy/NQvY6qSUI1NRJLSC1OCR3PUU2Mevx7MOYhBVwE5OIzzMJ+tRxGUsjSCvOqC0LQLVaXSipxjeJ9D46JwXi9+b0yh2XzUrbtzho3AkkMxrnDoUvC0KWHd46FIVH7/Iky6CYP5RsCXzezNwPzzezZwB8B/9bpwZ0E4ovN7CPAlwkdtgBw93/Z15KKDLud43mjtpqakQPzzCCO+U2sOWSplrd2cgo9i4u0camHNc3exbU8b/QyzrwYNxySvtUELHXMQzq4CMaJhfS10ewIFgKik+XFhCJxeFQYfMxodTQULi8tiRiXr6tZhdF9uCZGqN1mpT++W3dlHLIgBOPQyaz5x/nB7Y+SmlFNpg7xEhlQbwLOJ2SOX0Fok+7qMogvB34VqFKa5AdQIJaBcP/WnY0fzCOXL5p235n2tNVL+crGzcyrJFSThFFPqMaJNnKK2Z+DagL1emmsMCEI57Smo92bQ3ygNV0NHmvHOZknWB4m9ih6UycW2mSLIU/g1JtzdMSFIyCPQ50coJ4zkqaMpiMYEzhpyzSelX1ooUqstXc0xGk9rVF68nhDUV4Eox0tDjEchjQ1fRZwpbt/eH8O7iQQP9ndj9+fk4v02j0P7Wh0boLm1IprDz0IgPu2hKkZi05KMxG4t+2qsaCaMlpJgArzKwl50jp2smhLSi2kZrPM27YJ522amWqZ8/lPXcknP/A+fu/C1/HfN5xD5s5oWuoGnXjLVNEpRi13JjIaAbqeN9ubSVoD5ng9x1NjXjrSHMJUWhBi/NEdjC48aNrrsHNsV+iAlXuj1p3Eduxi+MqhCypsHqs3hlG1S33HCjMQ5vcuz2gms9PwxWHOBP7GzP4duAa43t07XqWwkwk/bzKz4/a3dCK90m6VosK9W3Y0gnDZT7fu7Ho5tu2qsX28zlgtZ8dEPaSBJ5UtjSnqIiA3JwEJwbiee9sgXPjkB97Hll9s4hMfeB+1LI/TYmaM13PquTNe9/i5zfmnMw9Bt0hpVxJjXiUhiXNlF+nwYv9aDuM55JaSWzqlDOM7t+/1WlQSY361aLu2OJVmeL59PAT4QxZMrQ+U2xCLIFxRiloGlLu/HDgW+CfChCE/ik26HekkED8D+K6Z3W1mt5vZHWZ2+/4VV6S7Htr+aJyVqvkHuxxQ2ikCz0x47j98g0fG64zVsvgIgbE8ZHdyb+DQX8kalcGiYphgYXaq2AGrmCIyTYyX/NFrOfgxh/PiV/5p6HiV5Y025LFaRi3z2AM7BPUi+Badv3bF/HTRbpuaMVoxqjFTkNO8IZjI8tB+bFaacSsu/vDIVsYf2dr2WixaMD98r/idl4ympOXOa9asGS2f1xroi2tUpKrL02MWC07ILFX6/TzQxyBx9xrwRUKN+FZCurojna6+JAI0l/4bhPTg5kfGgBCk3K0RPPK8NAWkx5X7aP4xb47BhVUzkJq+/WfbWb5whFXL57Ny0Si1RaNkk5pVi3RrUgqyFmvJTjNgZu6kMSAn5qRJSi1zznzJuTznRS9rBmGgmibNgJs6VU/IcqjG3k5FKrpQPJ9XSRhJjRxnJA3lqMYe3rnDwQe1LsDQLvDWNv+U6iFTFz8plo1cBIzveJilI6NYNkFeGeXh3Vnj5mn7eB5rzOGmIClq0MR/3RtrJMvsN2y9ps1sPWH6zN8mLJD0EeBFnR7fycxaP9nbPiK99ssdIQibxbmRS22RZUWv28SspfftTARggG1bw43K/JGUXROhVlzLQ0eqLHey0lCd4o9REstPEmqkWRziUwThLDSgkmJMZCGwXvfJK7ny7/+aF7/yTznld39/SjlCkM1ZUE3jdfDGzUCxglMRoOdXQ5q6miQtgXByAG7HPMcmxjq6NsXY4/FHdpHUdrNsdLQxfnjZvJRHJrJmEI69vosg3NLg3WapRpE++wNCTfgVcSWmfbLH1LSZ3ba3gzvZR6Sbtu4YY2sMwsUf7GJaxDQpAorFFCiNVGsvgjDAN998Co8+OsGuiaylXTjLw5Chemme56JcaWKNslcSYySN/8ZaaSWB0dQa7brVJOGqv/9rtjz4ANd88NLG8eWbkCLIFsG3CLwLRxIWjqTMryaMVoyDRhJG4meOpMa8+O90Q4Ysm2guj+g5PrIAH+l82cLRxcsZPWgp8+bPJ8Ebj3IQLr5TEYRtUkpcZq9wi9idx6Bw9w3u/q97CsJm9s3pjp+uRvz4vbQFG7CkgzLKkBmEtHSR1i2rxOkdRyu0pILLy+zNZBAufP11Jzeef+hbP2G8nrOg6kzkkMbhQqk3069FY9dIatSy0BZbpNUzL4Y0FasqhaD6B3/8ej76t+/lxa/808ZnFZNkZe5US9NgFCnoxaMpOd5on00wRtIw/rgSA18RgJNkz91HRpatZGLrA1jRi7oW5pmub7qHymFr9+laWW0XXhnFPKeapI3OXEBrOjoG4NFF+pMzDIYsM92JedO9OV0g/tUOTp7tfRcZJoMQhPc01jSxEMzqeagJT+7MccSyhTNfuDa2jE3wW8esAELv7dQMT5oLHRT/JmbkMf55nOUic6hZvLGIbd0jFeP5v3cOp539UnbHTlfl1WyKntDVJKGW58yrJMyrhBMnWKNXWyUJ16sIxBWDkU13gudMHL6XEYt5jnk9TIVZ29U6H/U+sGwirHWcZ8ynRq0SOnglTArCc/AvtwyVaX+A9xiI1TYsg6ioQLYLxkUteSQtOjvFIUG9LWKLVzz9qJbXj10Rxt5ufmSMzJup8wTHS7NuJUkSgmvjyzqJG5/82Ee59D1/xfmveQOnvPD3G/sXad1qrMlWSzXaYg1jaI4hNmAkTRop6SSrUd10V6PmObLpTnjsntdbL4Jksjv2YPacyuGP27+LlGeY52Sji6i4E9Z3LKWiJ7cRy6w3aD2e+63j9RJF+u3hnc1OQeVf5PJatkV8DrM3hYCcAocs7rwNsxeK8uzavbuRgg3DfIoAmmBxxaRa3hxj/H/eewm/2PQAl7//PTxvw8uAMPlGNU3iPs68JAwTyvIQoBeOJCRYYw7qYvhQNYHRBKwemrWSY07suPzVQ9dQe+jHeDoSarX7yerNY1vagyfVgssrPcns5qWb5Dlk2lHw+5dPEumTWh7aTSdyZyJ3xjNnV90bq/o0pmmMklgDfGj7o30r83Tmz5sX0v1ZDbIaVh8PQTmvhwAda8lFCvkNb3oTKw8/gle99o2MpikLquFR1IbDc4sTeXjLTF1GqIFXE2NBNWHRSEqyaxs28eh+tb1WD11D9TFHUzn8cftdG64eciQjy1aS7NpOuuMXpDseItm5BfKssfbx5JWeRAZNu0mvzOzk0suXTXd8J+sRvxq42t237XPpRLqoCCjFYgXlu+o0MYppMzJvnTSjvJgAwMol/Wkr7ojFuTD3cP/80nPP4/kvOTfchNQdSMhTb/TWGKtlLKimjdWaoLSGsRtViBOEGIlnkFQYWdrxsqkzpgjk9fvvxDwnmdhJtvgwACa2PsDI8sP7WTzpsiFsaLjWzK4C3k3omPVuYB3w6wDu/v3pDu6kRrwSuNnMrjWz9aZFQqUPHt45FhZFcG9JbYVJPCiNvQ2Pet6cUWpW/dIXadnieZQXtf3G2OPmTFxAnNc6DFsqZtfaXQ+zbRWLSBQrLxWrQJnneLX/ne/KKqueQLr6+PB4ZBMcQNpbBlcxpeqBPgbI04HVwDeAm4EHgN/o9OC9BmJ3/1/AWuBywqDle8zsL8zsmP0prci+KtqGi7jTXB4wBOC8tJZvWNig+ajHR3lBgwcHNE1dsHIwJqTaPaaYi9R7cUMCsbd1rAHvrueM1UIAHq/njTWHG+cqPZ+wyl4XbuindPXxjCw/XLVhmQ1qwC5gPqFGfJ975z0MO2oj9nAb/mB81IFlwKfN7N37XFyRA9Sux3QRhOuTAvGuujdqzAN2B91iSmck9zBrVWPN4abmXLveeNTiog4Qp7zMm9+7rDjVrlrOtt0afSi957T+7B7IY4DcTAjEv0ZYn+ElZvbpTg/upI34NcC5wBbC/JlvcPeahe6d9wB/tj+lFunE1h1jjV7RoUU0NIeahakiy6Gknjfbj8srF5klja7V6Z4GIQ+A0cXLm/M4ew5uodNWWiVNrBFoizbfMFtX3gjC442FHMIpUoup6swbY5UhzOs8qqWMpI8GKoR2x/nufkt8/iBwlpk1OmiZ2bLp+ll1MnxpBfA7k8cVu3tuZmfsT4lF9kURhBu9f62Ynzks9FCPy/3tKtpEG72FndE0JbVQG04TSAf0L8DEtgfBkpB+t8asHkCoFRe14KLmn+XOeNZsB67lTi1OJ5bGYU/FkCYI1+iXYzUOGg1zT689VD2RRbqlFITL264qvfwycMKeju9k0Ye3TfPeXXs7XmR/FXNKlxWLPDgh9Vosf5gTaoTZpJTV5PRV5t6YVGMQeVKZ8jrHyPI8zlPtjVWUisd4llPL8pZ28CQ1alnOaCVMDDJed5752IN7/G1E2puDE3pMm4LShB4y8Mq/s+4+ZT1foJGmLZ4XSwuW17HNclh76GAG4ZFlKxvPx3duB0twS/AkJWssFuHU8zB39O7YEauoCReTfhTfN3dn50TGYdVRxmpZY4rNYZT95HsApEc9uc8lkU4NVvNuT+zfFJcig6Q8ZaWXhiSFxRCaapk3glFqFgNyT4t6wEYXLWH7o7tIgSz2fp7ImsF4rBaWVxzPcnbXc+pZs214tBJ7UGfOhicf0d8vIiIdUSCWgVSkpYu24XLwhTgmNqamW6e7DMG3eJ5Y63jb2WLJwvk8uP3R0DErBuDddWciy9k5UQ+BuB7S0bU8J8+d0UrKeD3nZSes6nfxe0o14dklrCI256rESk3L7FTuoJUzNZ3l7nGSi/A6MSMJfZEateHUjNgFamDT0nsykXljopKJzNlVCxN17JjI2FXLGu3CRec0gP9x4lF7OatInw3eZBxdYWbPANa6+0fN7BBgkbvfF98+ZbpjFYil7zY9/GjL8BporQEXtd/ybFoeA1AxZGlXLWfhSNr4BQ9zNIdz1nLnuJWLe/RtDtzGzTtwp9Exa7yes2Oizng9BOKHHh1n+1iNiXrOTT/6JScdEzph/fnpnaxcKiLdZmYXE6a0fBzwUaAKfJw4u5a7b53ueAViGQj13BtL+jUD7fRLGRY31ROZx9WGQvtp0UZcTRIWjszOdU3KQbhIQ4e24ZyJek6WO1//wS9I04T/d/dmvvyaZ/a7yCIdG8Je0y8AngrcBuDuD5hZxyk4BWLpuywPPaGLwFv0im7XLlwohiXVcm/Ufo8/bN9XEBo09zwU1vdtDFHK80YP6fF6zng9I8ud+SMpnjv1POPrrzu5v4UW2QdhytZ+l6LrJtzdzcwBzGyfVpZRIJaB0LKIQ9bsymE051IuG8+KlZbCv8MUhIux0EUQLuaQLoYrpYnx2mcewxt+69g+l1hEomvN7EPAUjP7Q+A84MOdHqxALH2XTarpOmHMLxRzI0+9fS5SW8MQgO/dsqPxfXPCLGHFesIQav1hPukwYcnnv/MAr32m1lyR2WvYek27+3vN7NnAI4R24re5+w2dHq9ALH1175YdjX79zaFK4Rc1weIqQ+HfovPVbOv93KkiCOd4Yw7pRq/x+GSinlOvabEGmd2GLTVtZkcD/1EEXzObb2Zr3P3HnRw/O3uyyFC4+6FHgGJZv9YgnOXNpf5CjTDMJb2rPqtWF+5IMeVmljfXWS2GJBXzRyeJNVZTOmPd3BonLDIL/BOtfUqzuK0jqhFL34QUrLfMfFXUALPSGOFazF0/bfVwLVRQDFPK4qIO9UYK2huLORTCIg7GwQtHOG/dkX0stciBKZZBHDIVd58oXrj7hJmNdHqwasTSN8etXNwYplQ8oFkTdoexWkYtz4cuCBcdswrF36Wip/RkqcGS0Qo7xuu9KJ7IzHEa/SAO9DFANpvZmcULMzuLsHRwR1Qjlr6q5XnLwgxAYylDgJOOWt6PYvWEO+zO8pZafy3Pp/yBKZY0vG/b1NWoRGQgvBK42sz+jjDY42fAOZ0erEAsffMf9/6SeZWkJTVd1AaHOQBDaPdutIFn5c5Z3khVF8JEJ8YxyxdyxuMf068ii3TFMKam3f1HwElmtggwd9+xt2PKFIilZz52689IzFixoArAopEKY7EHcHms8DAv2Vc4buVibn9ge2O8cPHdi3HRWWlZw2pq/PyR8b6VVaS7mnOjz3Zm9lJ3/7iZvXbSdgDc/X2dnEeBWPpm2+4ao2nCeMzFnnncyr0cMXyKWnAt2/Mfp9SMI5fMmxM3KCKzTDGD1gGNqVQglp4592mrW15fddv9c27JvjIzGinpWu6NWjBAWsrXKwDLMBmm1LS7fyj+++fT7WdmF7n7X+7pffWalr6Zy0EYQmetoqNaYiH4VlMjTYzEwrZT1x7S51KKdNlw9prem7One1M1YpE+q6YWhv9TzB4W/lVNWGRo2HRvKhCL9ND9W3cCodd0WWLWXL4xnfZ3VmRWG6bU9D6Y9gsrEIv0yaKRpNFrPFUjkcwhw9Jreh9Me3etX3+RHsoJ82oXK0wdPL/KaJq27LNopMKiEd0jiwyRaeed1m+7SI/l7o3b4/JycEVqOjWlpmV4Faus9ZuZ/W/gD4HNcdOb3f0L+3muecD5wBOAecV2dz8v/vsX0x2vQCzSJ2kCKUaawKK00pjm8gmHLe5zyURmkNMyVK/PLnX393bhPFcBPwROB94O/D5wV6cHKzUt0mNFzdfiY8loSE2vWFBhxQLdG4vMQse6+1uBR939Y8DzgOM7PViBWKSHEsLc0dUkjBeuxsfyeSmpwVEHL+p3EUVmlBNmk+vGA1hhZreUHhfsY3FebWa3m9kVZrbsAL5WLf77sJk9EVgCrOn0YN1+i/TAT7funHLXW7wuZtE6bOlCROaCrHuZ6S3uvm5Pb5rZjUC7uXPfAnwQeAeh2fodwF8D5+1nOS6LgfytwHXAIuBtnR6sQCzSJ8XE8ArAIjPD3U/tZD8z+zDwuQP4nI/Ep18HHruvxys1LTLD7tuyA/cwdKmsmMZSZC4pJvToUmp6v5nZYaWXLwC+fwDneoyZXW5mX4yvjzOz8zs9XjVikRl075YdGHEmIYdqJaEYnbRyiWrCMgcNTq/pd5vZUwj3Bj8GXnEA5/pH4KOElDfAfwGfAi7v5GAFYpEuu2/LjkYbWGLht7yo+CoNLXPdoExx6e4v6+LpVrj7tWZ2UTx33cyyTg9WalpkBhU3/k5XO6iIyGB51MwOJs4pbWYnAds7PVg1YpEZkFoz8A5GFk5kcAzhTelrCb2ljzGz/wQOAV7Y6cEKxCJddN+WHY3nRfataBM+9pCD+lAikcEyKKnpbnL328zst4DHEVqi7nb32l4Oa1BqWqRLNm4uBWGgksBE5kxkriAsMsTM7GxgvrvfCTwf+JSZndDp8QrEIl1UTrntqjsLqsb8isYoiTS4k+fdeQyQt7r7DjN7BmG+6Y8RJgzpiFLTIl2UGNTzkI4eSY2jV6gmLFI2pB0Xix7SzwM+6O6fjas7daQvNWIze4+Z/TDO8fkZM1taeu8iM9toZneb2eml7U8zszvie++3OC2RmY2a2afi9m+Z2ZrSMeea2T3xcW5p+9Fx33visSNxu8Vzb4xl6zi1IOIeOmapTVhkzvm5mX0IeBHwBTMbZR/ia79S0zcAT3T3JxEGPl8EYTYSYANhTcf1wAfMrFg1/YPABcDa+Fgft58PbHP3Y4FLgUviuZYDFwNPB04ELi5N6n0JYfmrtcC2eA6A55TOfwH7kFoQyXHcQ0BWEBbZs0GYWavLXgRcD6x394eB5cAbOj24L4HY3b/k7vX48iZgVXx+FnCNu4+7+33ARuDEOBXZYnf/prs7cCWhQbw45mPx+aeBU2Jt+XTgBnff6u7bCMF/fXzvWXFf4rHlc13pwU3A0knToIns0c7xnF318BCR9kJq2rvyGBTuPubu/+Lu98TXm9z9S50ePwidtc4DvhifHwH8rPTe/XHbEfH55O0tx8Tgvh04eJpzHQw8XLoRaHuuNu+1MLMLiqW3Nm/e3NEXleE2v5pQSYwnHb6k30URkVlkxjprTbf8lLt/Nu7zFqAOXF0c1mZ/n2b7/hyzP+eautH9MuAygHXr1g3OrZn0xDd+/EsAqklCYkaahOUMjz9MQVhkWs6g9XjuuxkLxHtbfip2njoDOCWmmyHUQFeXdlsFPBC3r2qzvXzM/WZWISzIvDVuP3nSMV8DthBSzpVYK253rnafIwLAVzZuZl6lNZmU5fCUIxSERfZmSHtNH5B+9ZpeD7wRONPdx0pvXQdsiD2hjyZ0mvq2u28CdpjZSbGN9xzgs6Vjih7RLwS+EgP79cBpZrYsdtI6Dbg+vvdVmtOPnTvpXOfE3tMnAdvjZ4vwubt+wY33NJshqkn49cndedrqpXs6TERkWv0aR/x3wCiNdpZcAAAPxElEQVRwQxyFdJO7v9Ld7zSza4EfEFLWF7p7MT7rVYSlpuYT2pSLduXLgavMbCOhJrwBwN23mtk7gJvjfm93963x+RuBa8zsncB3aC5V9QXguYROYmPAy7v9xWV2+Nxdv2Dbrhq5O4kZ1dRYNNL+1+XXjlzWdruItDdgPZ77ri+BOA412tN77wLe1Wb7LcAT22zfDZy9h3NdAVzRZvu9hCFNk7c7cOF0ZZe5K409CHKH/7bm4P4WRmSWcgarx/MgGIRe0yID54zHP4aXnbCKn2wd475fPsqGJx/BWC1j++4ap649pN/FE5EhoikuRabxtmc/rvH8d48/vI8lERkSDpl6TbdQIBYRkZ5xFIgnU2paRESkj1QjFhGRnnGlpqdQIBYRkZ5SIG6l1LSIiEgfqUYsIiI947hqxJMoEIuISO+ojXgKpaZFRET6SDViERHpGY0jnkqBWEREekbDl6ZSalpERKSPVCMWEZGeUo24lQKxiIj0jIYvTaXUtIiISB+pRiwiIj3jDnXViFsoEIuISE8pNd1KqWkREZE+Uo1YRER6RuOIp1IgFhGRnspcgbhMqWkREZE+Uo1YRER6RuOIp1IgFhGRnlEb8VRKTYuIiPSRasQiItJTqhG3UiAWEZGeCesR5/0uxkBRalpERKSPVCMWEZHecfWankyBWEREeiakphWIy5SaFhER6SMFYhER6ZliGcRuPA6EmZ1tZneaWW5m6ya9d5GZbTSzu83s9AP6oA4oNS0iIj0zQKnp7wO/A3yovNHMjgM2AE8ADgduNLNfcfdspgqiGrGIiMw57n6Xu9/d5q2zgGvcfdzd7wM2AifOZFlUIxYRkd7p7hSXK8zsltLry9z9sgM85xHATaXX98dtM0aBWEREeqbLiz5scfd1e3rTzG4EVrZ56y3u/tk9HdZm24zm0hWIRURkKLn7qftx2P3A6tLrVcAD3SlRewrEIiLSUwPSWWtPrgM+YWbvI3TWWgt8eyY/UIFYRER6ZlCWQTSzFwB/CxwCfN7Mvuvup7v7nWZ2LfADoA5cOJM9pkGBWEREeswHIBC7+2eAz+zhvXcB7+pVWTR8SUREpI9UIxYRkZ5xh3wAasSDRIFYRER6yHFXIC5TalpERKSPVCMWEZGeGoTOWoNEgVhERHpHbcRTKDUtIiLSR6oRi4hIzzjgeb9LMVgUiEVEpKfUa7qVUtMiIiJ9pBqxiIj0jjprTaFALCIiPeQavjSJUtMiIiJ9pBqxiIj0TOg1rRpxmQKxiIj0jkOuXtMtlJoWERHpI9WIRUSkp5SabqVALCIiPaVA3EqpaRERkT5SjVhERHrG3TWhxyR9rRGb2evNzM1sRWnbRWa20czuNrPTS9ufZmZ3xPfeb2YWt4+a2afi9m+Z2ZrSMeea2T3xcW5p+9Fx33visSNxu8VzbzSz283shF5cBxGRucTdu/IYFn0LxGa2Gng28NPStuOADcATgPXAB8wsjW9/ELgAWBsf6+P284Ft7n4scClwSTzXcuBi4OnAicDFZrYsHnMJcKm7rwW2xXMAPKd0/gviZ4qIiMyYftaILwX+jDC+u3AWcI27j7v7fcBG4EQzOwxY7O7f9HAbdCXw/NIxH4vPPw2cEmvLpwM3uPtWd98G3ACsj+89K+5LPLZ8ris9uAlYGj9bRES6xPPuPIZFX9qIzexM4Ofu/r2YYS4cAdxUen1/3FaLzydvL475GYC7181sO3BwefukYw4GHnb3+nTnmvTepjbf4QJCrZkjjzxyr99ZRETAtejDFDMWiM3sRmBlm7feArwZOK3dYW22+TTb9+eY/TnX1I3ulwGXAaxbt04/VSIisl9mLBC7+6nttpvZ8cDRQFEbXgXcZmYnEmqgq0u7rwIeiNtXtdlO6Zj7zawCLAG2xu0nTzrma8AWQsq5EmvF7c7V7nNERKQLNI64Vc/biN39Dnc/1N3XuPsaQvA7wd0fBK4DNsSe0EcTOk192903ATvM7KTYxnsO8Nl4yuuAokf0C4GvxHbk64HTzGxZ7KR1GnB9fO+rcV/iseVznRN7T58EbI+fLSIi3eAhEHfjMSwGahyxu99pZtcCPwDqwIXunsW3XwX8IzAf+GJ8AFwOXGVmGwk14Q3xXFvN7B3AzXG/t7v71vj8jcA1ZvZO4DvxHABfAJ5L6CQ2Brx8Jr6niIhIoe+BONaKy6/fBbyrzX63AE9ss303cPYezn0FcEWb7fcShjRN3u7AhR0WXURE9plr9aVJ+h6IRURk7tB6xFNprmkREZE+Uo1YRER6x1UjnkyBWEREekoTerRSalpERKSPVCMWEZGeGqaVk7pBgVhERHrGfbgm4+gGpaZFRET6SDViERHpKXXWaqVALCIiPeV5tved5hClpkVERPpINWIREekdd9WIJ1Eg7oJbb711i5n9pAunWkFYL1madE2m0jVpT9dlqm5dk6O6cA4AHAXiyRSIu8DdD+nGeczsFndf141zDQtdk6l0TdrTdZlK12R2UCAWEZHecfBMNeIyBWIREekhpaYnU6/pwXJZvwswgHRNptI1aU/XZSpdkz0ws7PN7E4zy81sXWn7GjPbZWbfjY9/mOmyqEY8QNxdvzST6JpMpWvSnq7LVAN5TQan1/T3gd8BPtTmvR+5+1N6VRAFYhER6alBCMTufheAmfW7KEpNi4iITHK0mX3HzL5uZs+c6Q9TIO4iM3u9mbmZrShtu8jMNprZ3WZ2emn708zsjvje+y3elpnZqJl9Km7/lpmtKR1zrpndEx/nlrYfHfe9Jx47ErdbPPdGM7vdzE7oxXWIn/0eM/th/NzPmNnS0ntz8pp0g5mtj9dto5m9qd/l2R9mttrMvmpmd8U2uj+J25eb2Q3x/+wGM1tWOqZvPzO9ZGZpDACfi6+H7poU44i78QBWmNktpccFk67njWb2/TaPs6Yp4ibgSHd/KvBa4BNmtvhAv/d0FIi7xMxWA88GflradhywAXgCsB74gJml8e0PAhcAa+Njfdx+PrDN3Y8FLgUuiedaDlwMPB04Ebi49Et5CXCpu68FtsVzADyndP4L4mf2yg3AE939ScB/ARfBnL8mByRep78nfIfjgJfE6znb1IHXufvjgZOAC+P3eBPw5fh/9uX4ehB+ZnrpT4C7Sq+H75o43QzEW9x9XenR0ibu7qe6+xPbPD67x+K5j7v7L+PzW4EfAb9ywN97GgrE3XMp8GdAeVmRs4Br4n/sfcBG4EQzOwxY7O7f9LBC9pXA80vHfCw+/zRwSryjPR24wd23uvs2QqBbH997VtyXeGz5XFd6cBOwNH72jHP3L7l7Pb68CVhVKtOcvCZdcCKw0d3vdfcJ4BrC95lV3H2Tu98Wn+8gBJ4jaP1/nvx/1s+fmZ4ws1XA84CPlDbP6WvSD2Z2SHFTY2aPJdzM3DuTn6lA3AVmdibwc3f/3qS3jgB+Vnp9f9x2RHw+eXvLMTGQbQcOnuZcBwMPl4Je23O1ea+XzgO+GJ/rmuy/2Vz2tmJ69KnAt4DHuPsmCMEaODTu1u+fmV75G8LNfF7aNoTXxMnzrCuPA2FmLzCz+4FfBz5vZtfHt34TuN3Mvke4CXmlu289oA/bC/Wa7pCZ3QisbPPWW4A3A6e1O6zNNp9m+/4csz/n6orprkmR+jGztxBSkVfvpUxDcU1m2Gwu+xRmtgj4Z+B/uvsjtufeq/3+mZlxZnYG8JC732pmJ3dySJtts+OaDMjwJXf/DPCZNtv/mfBz2TMKxB1y91PbbTez44Gjge/FPySrgNvM7ETCHeTq0u6rgAfi9lVttlM65n4zqwBLgK1x+8mTjvkaYUL3pWZWiXeu7c7V7nMO2J6uSSF2+jgDOCWmyaYr01Bckxk2m8vewsyqhD92V7v7v8TNvzCzw9x9U0yxPhS39/tnphd+AzjTzJ4LzAMWm9nHGcJr4gzG8KVBotT0AXL3O9z9UHdf4+5rCD/UJ7j7g8B1wIbYW/FoQlvDt2OKaYeZnRTbYc4Bis4D1wFFr8UXAl+JQex64DQzWxY7UpwGXB/f+2rcl3hs+VznWHASsL1Ic800M1sPvBE4093HSm/N2WvSBTcDa2NP1hFCZ53r+lymfRb/fy8H7nL395XeKv8/T/4/6+fPzIxz94vcfVX8G7IhlvelzOFrMpeoRjyD3P1OM7sW+AEhPXuhuxe3gq8C/hGYT2g/LdpQLweuMrONhLvVDfFcW83sHYQ/xgBvL7VbvBG4xszeCXwnngPgC8BzCR05xoCXz8T33IO/A0aBG2Km4CZ3f+UcvyYHxN3rZvZqwh/OFLjC3e/sc7H2x28ALwPuMLPvxm1vBv4KuNbMzieMPjgbBuL3qJ+G75q4a9GHSayZMRQREZlZlSVH+OJff0VXzrXt+otv9SFY5lGpaRERkT5SalpERHpnQHpNDxIFYhER6SEF4smUmhYREekj1YhFRKRnwjjifK/7zSUKxCIi0jtqI55CqWmRHjKzNWa2qzR+9kDP95Q4G1NXWFiecKeZzfohISKzhWrEIr33I3d/SpfO9RRgHWGiko6Upiycwt1/28y+1qWyibSlGnEr1YhFusTMfs3MbjezeWa20MKi90/cyzFrzOyHZvYRCwuWX21mp5rZf1pYjP3EuN9CM7vCzG62sHD8WXGay7cDLzaz75rZi9vtF4//AzP7JzP7N+BLZnaYmf17PO77ZvbMGb9AIgA+GKsvDRLViEW6xN1vNrPrgHcSphf8uLt/v4NDjyVMXXgBYZrB3wOeAZxJmPrx+YRVvr7i7ueZ2VLg28CNwNuAde7+agAz+4vJ+1lYJQvCcm9PilMavo4wn/C7LKy9uqAb10BE9p0CsUh3vZ0QTHcDr+nwmPvc/Q4AM7sT+LK7u5ndAayJ+5xGWJ3n9fH1PODINueabr8bSnMI3wxcYWEVpH919660WYvsjYPmmp5EgViku5YDi4AqIQg+2sEx46Xneel1TvN31IDfdfe7ywea2dMnnWu6/Rplcfd/N7PfBJ5HWAjgPe5+ZQdlFTkw6jU9hdqIRbrrMuCtwNXAJV087/XAH8el7TCzp8btO4CDOtivhZkdRViI/sOEFXVO6GJZRWQfqEYs0iVmdg5Qd/dPxHbXb5jZs9z9K104/TuAvwFuj0H2x8AZhPVi3xSHQ/3lNPtNdjLwBjOrATsJ69aK9IBqxJNpGUSRHjKzNcDn3H3a3tT9FIcvvd7db+l3WWT4JAsP8eqvPr8r55q47SNaBlFE9lkGLOnWhB7dZmZfBR4L1PpdFpG5QqlpkR5y958Bq/tdjj1x99/udxlkyKmz1hRKTYuISM+Y2f8FVnTpdFvcfX2XztU3CsQiIiJ9pDZiERGRPlIgFhER6SMFYhERkT5SIBYREekjBWIREZE++v92J66NyPvhyQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "ds.sel(time='2018-11-05').u.plot(ax=ax)\n",
    "ax.scatter(x, y, c='k', s=3)\n",
    "ax.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we do a nearest neighbour interpolation to get the closest u and v to x and y.\n",
    "\n",
    "To do this x and y need to be xarray DataArrays with the **same** dimensions.  The dimension name used for `ind_x` and `ind_y` does not have a dimension in the DataSet you are trying to subset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[85.039    85.916435 87.03461 ]\n",
      "[0.       6.340192 8.746162]\n",
      "[[-2.  -4.4 -4.9]]\n",
      "[[-6.1 -5.6 -4.2]]\n"
     ]
    }
   ],
   "source": [
    "ind_x = xr.DataArray(x, dims=['i'])\n",
    "ind_y = xr.DataArray(y, dims=['i'])\n",
    "print ( ds.sel(time='2018-11-05').latitude.sel(x=ind_x, y=ind_y, method='nearest').values )\n",
    "print ( ds.sel(time='2018-11-05').longitude.sel(x=ind_x, y=ind_y, method='nearest').values )\n",
    "print ( ds.sel(time='2018-11-05').u.sel(x=ind_x, y=ind_y, method='nearest').values )\n",
    "print ( ds.sel(time='2018-11-05').v.sel(x=ind_x, y=ind_y, method='nearest').values )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interoplate ice motion vectors to arbitrary points\n",
    "\n",
    "I use `scipy.interpolate.interp2d` for this.  This function returns 2d arrays.  So I use list comprehension to find individual x, y pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = ds.sel(time='2018-11-05').u.values\n",
    "u[np.isnan(u)] = -99.9\n",
    "v = ds.sel(time='2018-11-05').v.values\n",
    "v[np.isnan(v)] = -99.9\n",
    "fu = interp2d(ds.x, ds.y, u, kind='linear')\n",
    "fv = interp2d(ds.x, ds.y, v, kind='linear')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u =  [-1.94814917 -4.22065543 -4.93430168]\n",
      "v =  [-6.13456718 -5.56727877 -4.18583824]\n"
     ]
    }
   ],
   "source": [
    "print ('u = ', np.diag(fu(x, y)))\n",
    "print ('v = ', np.diag(fv(x, y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now write this as a nice function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def interp_uv(x, y, u, v, xi, yi):\n",
    "    \n",
    "    # Set NaN to -99.99 in u and v arrays\n",
    "    # - probably not the best approach but it will work for now\n",
    "    um = np.where(~np.isnan(u), u, -99.9)\n",
    "    vm = np.where(~np.isnan(v), v, -99.9)\n",
    "    \n",
    "    fu = interp2d(x, y, um, kind='linear')\n",
    "    fv = interp2d(x, y, vm, kind='linear')\n",
    "    \n",
    "    ui = np.array([fu(x_i, y_i) for x_i, y_i in zip(xi, yi)])\n",
    "    vi = np.array([fv(x_i, y_i) for x_i, y_i in zip(xi, yi)])\n",
    "    \n",
    "    return ui.flatten(), vi.flatten()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...and we get the same results yea!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.94814917 -4.22065543 -4.93430168]\n",
      "[-6.13456718 -5.56727877 -4.18583824]\n"
     ]
    }
   ],
   "source": [
    "ui, vi = interp_uv(ds.x, ds.y, \n",
    "                   ds.sel(time='2018-11-05').u.values, \n",
    "                   ds.sel(time='2018-11-05').v.values, \n",
    "                   x, y)\n",
    "print (ui)\n",
    "print (vi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
